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BULK  PROPAGATION  CHARACTERISTICS  OF  DISCRETE  RANDOM  MEDIA 

V.  N.  Bcingi,  T.  A.  Seliga,  V.  K.  Varadan  and  V.  V.  Varadan 
Wave  Propagation  Group,  Boyd  Laboratory 
The  Ohio  State  University,  Columbus ,  Ohio  43210 


The  propagation  of  electromagnetic  waves  in  an  infinite  med¬ 
ium  composed  of  a  random  distribution  of  identical,  finite 
scatterers  is  studied.  The  T-matrix  of  a  single  isolated 
scatterer,  obtained  by  using  the  null  field  equation,  is 
used  to  make  the  equations  for  the  field  incident  on  a  par¬ 
ticular  scatterer  and  the  field  scattered  by  it,  self 
consistent.  The  method  we  propose  is  well  suited  for  com¬ 
putations  at  wavelengths  comparable  to  obstacle  size  and  for 
non-spherical  obstacles.  The  attenuation  associated  with 
the  coherent  field  as  predicted  by  our  computations  is 
compared  with  the  only  two  sets  of  experimental  results 
that  can  be  found  in  the  literature.  Agreements  and  dis¬ 
crepancies  are  examined  and  the  range  of  validity  of  the 
assumed  quasi-crystalline  approximation  ISAM  is  discussed. 

Further  improvements  using  the  coherent  potential  approxi¬ 
mation  (CPA)  and  the  'self  consistent  approximation’  MCA, 
as  well  as  improved  models  of  the  pair  correlation  function 
are  suggested. 

A 

1.  INTRODUCTION 

A  treatment  of  the  propagation  of  electromagnetic  waves  in  an 
infinite  medium  composed  of  a  random  distribution  of  identical  finite 
scatterers  is  presented.  The  coherent  or  average  wave  in  such  a  med¬ 
ium  will  be  assumed  to  be  a  plane  wave  propagating  in  a  homogeneous 
continuum  characterized  by  a  "bulk"  complex  wave  number.  This  wave 
number  will  depend  on  both  the  frequency  and  the  concentration  of 
the  discrete  scatterers  causing  the  effective  medium  to  be  dispersive. 
The  aim  of  this  work  is  to  present  a  multiple  scattering  theory  and 
a  computational  method  of  obtaining  such  dispersion  relations  for 
random  media  models  including  comparisons  with  seme  past  laboratory 
experiments. 

The  effects  of  multiple  scattering  on  the  coherent  wave  are  of 
great  practical  importance ,  in  particular  the  dependence  on  concen¬ 
tration  at  wavelengths  comparable  to  scatterer  size.  At  very  low 
concentrations  (<  II  by  volume)  multiple  scattering  can  be  neglected 
and  each  scatterer  can  be  treated  as  independent.  However,  m  many 
practical  situations  the  concentration  can  range  between  11  to  201 
where  multiple  scattering  effects  may  be  important.  This  is  parti¬ 
cularly  reflected  in  the  study  of  higher  order  statistics  of  the  ran¬ 
dom  medium  including  radiative  transfer  theory  which  assume  that  the 
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coherent  wave  propagates  with  a  wave  number  that  is  only  dependent 
on  the  forward  scattering  amplitude  of  a  single  scatterer.  This 
result  is  obtained  as  a  solution  of  the  Poldy-Twersky  integral  equa¬ 
tion  for  an  infinite  slab  medium  filled  with  large  tenuous  scatter- 
ers  ill.  However,  this  well  known  analysis  neglects  correlation 
oetween  the  scatterers,  and  recent  advances  in  multiple  scattering 
theory  by  Twerskv  [2,3,41  have  accounted  for  the  effects  of  pair- 
correlation  at  higher  concentrations. 

The  theoretical/computational  method  presented  here  is  based  on 
a  self-consistent  exciting  field  approach  and  relies  on  the  T-matrix 
[5i  which  characterizes  the  response  of  an  isolated  single  scatterer 
to  an  arbitrary  exciting  field.  The  random  medium  is  described  sta¬ 
tistically  with  respect  to  the  random  positions  that  each  scatterer 
can  occupy  through  the  first  and  second  order  probability  distribu¬ 
tion  functions.  Ensemble  or  configurational  averaging  together  with 
Lax's  quasi-crystalline  approximation  !QCA)  yields  a  set  of  “hole- 
correction  integrals.  3y  assuming  a  plane  wave  behavior  for  the 
coherent  wave  and  using  the  extinction  theorem  gives  rise  to  a  homo¬ 
geneous  system  of  equations  whose  singular  solutions  yield  the  com¬ 
plex  wave  number.  The  method  is  necessarily  computational;  however, 
analytical  forms  of  the  dispersion  relations  are  obtained  in  the  low 
frequency  or  Rayleigh  limit  for  spherical  and  spheroidal  scatterer 
geometries.  This  paper  closely  follows  the  developments  given  by 
the  authors  [6,7,81  for  acoustic,  electromagnetic  and  elastic  waves. 
Previous  work  which  forms  the  basis  for  the  present  analysis  is 
given  in  Refs.  [9-12],  The  numerical  results  obtained  for  a  random 
distribution  of  spheres  are  compared  with  the  experimental  results 
of  Hawley,  Beard  and  Twersky  [13]  and  Olsen  and  Kharadly  [141. 

One  of  the  aims  of  this  paper  is  to  suggest  improvements  to  the 
QCA  by  incorporating  the  coherent  potential  approximation  (CPA) . 
Certain  multiple  scattering  processes  that  are  neglected  by  the 
OCA-CPA  scheme  can  be  restored  by  making  the  'self-consister.t  approx¬ 
imation'  (SCA) .  If  these  improvements  are  used  with  more  realistic 
forma  of  the  pair  correlation  function,  then  it  may  be  possible  to 
extend  the  present  formalism  to  a  wider  range  of  frequencies  and 
concentrations. 


2.  MULTIPLE  SCATTERING  FORMALISM 

Consider  N  identical  scatterers  located  in  free  space  as  shown 
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in  Fig.  1  wuere  0^  and  0,  refer  so  the  centers  of  the  i-th  and  j-th 
scaeterers.  The  scatterers  are  assumed  to  be  bodies  of  revolution 
with  parallel  symmetry  axes,  and  the  coherent  wave  is  assume  to 
propagate  along  this  direction  so  that  the  bulk  medium  is  isotropic 
and  does  noc  cause  any  depolarization.  Let  be  the  relative  die¬ 
lectric  constant  of  the  homogeneous  scatterers. 


Me  represent  an  incident  plane  electromagnetic  wave  propagating 
along  the  positive  z  axis  wich  wave  number  k  and  of  unit  amplitude 
with  an  e~li“  time  dependence  by 

E°(r)  -  e  exp(ikz)  tli 

where  e  is  the  unit  polarization  vector.  With  no  loss  of  generality 
we  choose  e  to  be  either  x  or  y. 

The  total  electric  field  at  any  point  in  free  space  outside  the 
scatterers  is  the  sum  of  the  incident  field  and  the  fields  scattered 
by  all  the  scatterers.  This  is  written  as 

-o  -  N  . 

E(r)  »  £°  !r)  *  '  Es(r  -  r.  )  .  12! 

i=l 

where  E*(r  -  rj  is  the  field,  scattered  by  the  i-th  scatterer,  at 
the  point  of  observation  r.  The  field  that  excites  the  i-th  scat¬ 
terer,  however,  is  the  incident  field  plus  the  fields  scattered  from 
all  the  other  scatterers.  The  exciting  field  term  Ee  i3  used  to 
distinguish  between  the  field  actually  incident  on  a  scatterer  and 
the  external  incident  field,  E°,  produced  by  a  source  at  infinity. 
Thus,  at  a  point  r  in  the  vicinity  of  the  i-th  scatterer,  we  write 

E®(r)  -  i°(r)  +  "  E?(r  -  r.)  ,  a  <  ,r  -  r.  !  <  2a  (3) 

1  if>i  3  3  1 

where  'a'  is  the  radius  of  the  imaginary  sphere  circumscribing  a 
scatterer  (see  Fig.  1).  In  this  analysis,  we  have  assumed  that 
there  is  no  interpenetration  of  the  imaginary  spherical  shells  of 
radius  'a'  which  circumscribe  each  scatterer. 

The  T-matrix  formulation  of  scattering  we  adopt  here  is  based  on 

the  extended  integral  equation  approach  due  to  Waterman  [151.  The 

scattered  and  exciting  fields  with  respect  to  a  particular  scatterer 

are  expanded  in  terms  of  a  complete  3et  of  basis  functions 

(M _ ,M  _ )  which  are  the  vector  spherical  harmonics.  These  form 

visn  t.  inn 
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the  total  field  at  any  point  in  free  space,  the  expansion  coeffi¬ 
cients  of  the  field  scattered  by  the  ]-th  scatterer  may  be  formally 
related  to  the  coefficients  of  the  field  exciting  the  j-th  scatterer 
through  the  T-matrix: 
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where  summation  is  denoted  by  the  repeated  index  convention.  The 
elements  of  the  T-matrix  involve  surface  integrals,  which  can  be 
evaluated  in  closed  form  for  spherical  geometry,  while  for  scatterers 
of  arbitrary  shape  they  can  only  be  evaluated  numerically.  The 
T-matrix  for  a  single  scatterer  is  of  the  form 


T  »  (Q-1)  ReQ 


18) 


where  ReQ  and  Q  are  matrices  which  are  functions  of  the  surface  S  of 
the  scatterer  and  of  the  nature  of  the  boundary  conditions. 

Substitution  of  Eqs.  (5)  and  (6)  in  Eq.  (3)  yields 


*  p  , 

"  '  7  'b*  ReM  (r-r.)  +  c1  ReN  (r  -  r, )  r  •  E°(r) 

pio  mio  _“e  '  -mp  ‘mp  1  !np  -“P  1  / 
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Note  that  the  series  on  the  right-hand  side  of  Eq.  (9)  is  expressed 
with  respect  to  the  center  of  the  j-th  scatterer  while  the  series  to 
the  left  is  expressed  with  respect  to  the  center  of  the  i-th  scat¬ 
terer.  The  addition  theorem  for  the  vector  basis  functions  will  be 
used  to  express  the  right  hand  side  of  Eq.  (9)  in  terms  of  the  center 
of  the  i-th  scatterer.  Formulas  for  the  translation  of  (MCr-r^), 
N(r-r^);  to  an  origin  centered  around  the  i-th  scatterer  can  be 
found  in  [16,171: 
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where  the  geometry  is  shown  in  Fig.  2.  )Jote  that  the  minimum  value 
of  r.^  is  2a  and  ' r-r^!  <  -a  ‘ or  Eq.  (6>  to  be  valid  so  that  the 


condition  for  use  of  Eq.  (IQ)  is  always  satisfied. 


=0 , 


It  then  remains  to  expand  the  incident  field  E  (r)  in  terms  of 
an  origin  centered  at  the  i-th  scatterer: 


E°(r)  -  x  expdk;^  exp  [ ikz •  ( r-rJ 


OQ  , 

!  ^  sh  \fol«  ^0l.««i>  +  *els  ReNels  (?-?,»  ;  Ul) 


where  I ^  -  ri'Z;  fols'  9els  are  the  '<nown  expansion  coefficients 

and  e  *  x.  Substituting  Eqs .  (10a),  (10b)  and  (11)  into  Eq.  (9) 
gives 
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Equation  .12!  is  but  an  expansion  of  Eq.  >3)  in  terms  of  vector 
spherical  harmonics  with  respect  tc  an  origin  3,  centered  at  the 
i-th  scacrerar.  A  relation  between  the  unxnown  expansion  coeffi¬ 
cients  ;b*  ,cl  )  of  the  exciting  field  and  the  unknown  expansion 
~mp  -p.p 

coefficients  iB".  ,C^  !  of  the  scattered  field  can  be  obtained 

j  •  n  j  »-  n 

from  Eq.  *12)  by  using  the  known  orthogonality  properties  of  the 
vector  spherical  harmonics.  It  may  be  shown  that 
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where  [  denotes  the  sum  over  all  scatterers  except  the  i-th 
scatterer . 


Substituting  Eq.  (7)  into  the  above  expressions,  we  obtain  a 
self-consistent  set  of  equations  for  the  unknown  exciting  field 
coefficients  as  given  in  Eqs.  !14a,b). 
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Thus,  the  unknown  scattered  field  expansion  coefficients  are  elimin¬ 
ated  through  use  of  the  T-matnx  resulting  m  a  set  of  equations 
mvolvina  the  expansion  coefficients  ib  , c  )  of  the  excitmo 

-.Tip 

field  only.  These  coefficients  are  functions  of  the  positions  of 
all  the  scatterers. 


3.  CONFIGURATIONAL  AVERAGING 

In  order  to  average  the  wave  fields  over  the  positions  of  all 
the  scatterers  we  define  the  probability  density  function  of  finding 
the  first  scatterer  at  r,  ,  the  second  scatterer  at  r,  ,  and  sc 
forth  by  p i r^ , r ^ >  • • • / ) .  This  probability  density  function  may  oe 
written  as 

p(r:,r:,  •  •  •  ,r^)  -p(rj  ?  [  r: ,  r2  r,. 1  rj 

»  p(r.)  p ir,r, ) (p (r, , r^ 1 rv  r,,r.) 


where  ?(r.)  denotes  the  probability  density  of  finding  a  scatterer 
at  r^  ,  while  p(r^  r^  denotes  the  conditional  probability  of 
finding  a  scatterer  at  r^  if  a  scatterer  is  <nown  to  be  at  r. .  A 
prime  in  the  first  expansion  of  Eq.  (13)  means  is  absent  while 

two  primes  in  the  second  expansion  of  Eq.  (13)  means  both  r.  and 
r^  are  absent. 

If  the  scatterers  are  randomly  distributed,  the  positions  of  ail 
scatterers  are  equally  probable  within  the  volume  V  accessible  tc 
the  scatterers,  and  hence 


P<rA> 


where  n  is  the  uniform  number 
o 

the  total  number  of  scatterers. 
imaginary  spheres  circumscribing 
conditional  density  as  follows: 


density  of  the  scatterers  ar.d  *J  is 
In  addition,  for  nonoverlap  of  the 
each  scatterer  we  approximate  the 
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The  f  o  r  m  of  the  pair '“correlation  in  Zq .  ,17)  describes  the  usuai 
radially  symmetric  distribution  function  with  an  exclusion  surface 
or  "hole"  corresponding  to  a  sphere  of  radius  2a.  This  form  of  tr.e 
pair-correlation  is  expected  to  be  satisfactory  at  low  concentra¬ 
tions  c  where  c  =  n  ^  y  ~a  "  [13].  However,  at  hi  crier  concentra¬ 
tions  Eq.  il?)  dees  not  account  for  either  the  space  occupying 
property  of  the  sea t terers  no r  the  increase  m  local  order  wit n 
increasing  c  resulting  from  the  well  defined  shape  of  tr.e  neigh¬ 
boring  scatterers  [21 .  The  available  volume  V  for  locating  the 
center  of  the  S-th  scatterer  randomly  m  the  total  volume  V  after 
tr.e  first  N-i  scatterers  have  been  previously  located  is  a  func¬ 
tion  of  concentra t i on .  For  hard  spheres,  the  minimum  value  of 
’’  =  V  -  N’-l'  ~'2a,|J  =  V  i  l-3c)  .  The  maximum  value  is  consideracl 

higher  than  this.  Since  this  restriction  applies  only  to  the  cen¬ 
ters  of  spheres  and  dees  not  taxe  into  account  their  physical  dimen¬ 
sions  ,  this  model  is  expected  to  be  better  in  the  Rayleigh  regime 
when  the  spheres  appear  to  be  point  particles  :see  Section  **)  .  Thus 
a  realistic  form  of  the  pair  correlation  function  must  depend  on 
concentration  as  well  as  the  distance  between  the  two  scatterers 
under  consideration.  Twer  sky  [2]  r.as  considered  concentration 
dependent  pair  correlations  and  has  extended  the  formalism  to  higher 
concentrations  in  the  Rayleigh  regime. 

ve  denote  the  configurational  average  of  a  statistical  quantity 
5  as 
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where  m  Zc.  (19a)  we  have  averaged  over  all  scatterers  except  the 
l-th  and  in  Eq.  vl9b),  over  all  scatterers  except  the  i-th  and  i-th, 
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and  so  on.  Multiplying  noth  sides  of  Ecs  .  (14a, b)  by  the  probabili¬ 
ty  density  given  by  Eq.  (13)  and  using  Eqs.  '  1 6  )  —  ■•  1 3  )  ,  we  obtain  tr.e 
configurational  average  of  b  ^ 
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where  V  denotes  the  volume  of  the  medium  excluding  the  volume  of 

N  ( 

tne  nole  of  radius  2a.  For  identical  scatterers,  ”  can  be 

j  *  1 

replaced  by  CM— 1 ) .  Equations  (19,20)  indicate  that  the  conditional 
average  with  one  scatterer  fixed,  vk.  ,  { •- t>irnp> i  <  ‘•c*  ls 

given  in  terms  of  the  conditional  average  with  two  scatterers  fixed, 

viz.,  [cb-'  »  ■  ,  ‘C1  '■■].  Lax  [10)  has  suggested  a  cuasi- 
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crystalline  approximation  (CCA)  to  close  the  system: 
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The  validity  of  this  approximation 


j 


is  examined  in  Section  3. 
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4 .  THE  COHERENT  WAVE 


The  heirarchy  of  equations  implied  by  Ec3 .  (19,20)  is  truncated 
by  invoking  the  CCA,  i.e.,  Eq.  (21).  A  plane  wave  solution  for  the 
system  of  equations  given  by  Eqs.  .19)  and  <20)  is  assumed,  using  ar. 
effective  wave  number  X  to  characterize  the  bulk  medium : 


,i 

w  ..sip  i 


.  p 


k  iK  •  r 


-mp 


i 


are  unknown  constants.  The  effective  wave 
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number  K  is  assumed  to  be  parallel  to  that  of  the  incident  wave 
w.-.irn  m  the  present  case  is  along  the  z-axis.  Since  the  symmetry 
axes  of  the  scatterars  are  also  assumed  parallel  to  the  z-axis,  the 
ou-<  medium  is  isotropic.  Substitution  of  Eq.  ;22i  in  Eqs.  ,19'  and 
.21  gives 
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where  i.  is  the  Kronecker  delta  and  applies  to  the  T-matrix  ele- 

«.r. 

ments  of  rotationally  symmetric  scatterers,  viz.,  no  azimuthal  mode 
coupling . 
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It  remains  to  perform  trie  integration  over  V  ir.  Ecs.  23a, c;, 
the  details  of  which  are  given  m  Appendix  A.  To  derive  tne  disper¬ 
sion  relations,  we  apply  the  extinction  theorem  to  the  two  sets  cf 
terms  m  Eq.  i23)  after  integration  (each  satisfying  the  wave  equa¬ 
tion  with  wave  numbers  K  and  k)  as  discussed  in  detail  by  Tver sky 
[4]  and  Varadan,  Varadan  and  Pao  [31.  According  to  the  extinction 
theorem  the  integral  1^  which  is  evaluated  over  a  surface  So 
precisely  cancels  the  incident  field  (see  Appendix  A!.  3v  equating 
the  remaining  terms,  we  obtain  an  infinite  system  of  equations  for 
the  unknowns: 
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We  define  c  as  the  effective  "spherical"  concentration  which 
equals  4-a3no/3.  The  term  (JH) x  is  given  by 

(JH>,  «  2ka  jx(2Ka)  h^(2Ka)  -  2Ka  h. (2ka)  j‘(2Ka)  .  125) 

The  factors  *  _(n<n',l)  and  x  (n,n',U  are  oiven  in  Appen- 
oo  oe 

dix  3.  They  are  related  to  those  given  by  Cretan  [17J  which  are 
also  summarized  in  Appendix  3. 

The  set  of  equations  given  in  Eq.  (24a, b)  are  homogeneous  and 

linear  in  the  unknowns  (Y  ._,Z  .  i.  Tor  a  nontrivial  solution,  we 

oin  eln 

require  that  the  determinant  of  the  coefficient  matrix  vanish  which 
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yields  a  relation  for  the  effective  wave  number  X  m  terms  of  k 
and  the  T-matrix  of  the  scatterer.  This  is  the  dispersion  relation 
for  the  scatterer-f llled  medium.  Equations  (24a)  and  (24b)  form  a 
general  expression  valid  for  any  arbitrary  collection  of  scatterers 
orovided  the  scatterers  are  identical  and  rotationally  symmetric 
with  parallel  orientation  along  the  X  vector.  Since  the  T-matrix 
is  the  only  factor  that  contains  information  about  the  exact  shape 
and  boundary  conditions  at  the  scatterers,  one  can  also  use  the 
above  formalism  for  a  collection  of  perfectly  conducting,  dielectric 
or  multi-layered  scatterers.  The  T-matrix  for  such  various  scat¬ 
terers  nas  been  studied  by  many  authors  [23,24,25]. 


5.  LOW  FREQUENCY  SOLUTIONS 

In  the  Rayleigh  or  low-frequency  limit,  the  sire  cf  the  scat¬ 
ters  -  is  considered  to  be  very  small  compared  to  the  incident  wave- 
lengu...  It  is  then  sufficient  to  take  only  the  lowest  order  coef¬ 
ficients  m  the  expansion  of  the  fields.  In  this  limit,  the  elements 
of  the  T-raatnx  can  be  obtained  in  closed  form  for  simple  shapes 
such  as  sphere  and  spheroid  [221: 
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where  e  is  defined  by  e  »  .  (a  'b)  -1  for  the  oblate  spheroid  i’a1 
and  'b'  being  the  semi-major  and  semi-mir.or  axes,  respectively). 

The  functions  f,  and  fj  are  given  by 


f  2  (e ) 


'23a) 


From  Eq.  (24b)  and  using  leading  terms  of  the  T-matrix  of  0(k3a3), 
we  obtain  the  following  result  for  the  unknown  Zell! 
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Dispersion  relations  ire  obtained  by  substituting  Eq.  (26)  or  (2 
m  Eq.  i. 2 9 1  and  using  the  leading  tern  m  the  expansion  for  the 
Sessei  and  Hankel  functions  composing  ( JH>,  and  ( JH),: 
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Equation  (30)  is  recognized  as  the  dispersion  relation  of  the 
Clausius-Mossotti  form. 


If  the  concentration  c  <<  1  ,  the  dispersion  relations  simplif 
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Equation  (33)  can  be  written  ir.  terms  of  the  forward  scattering  ampli 
tude  F  ( 0 )  : 


5? 
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which  in  this  case  is  valid  only  for  very  low  concentrations 
tne  Rayleigh  limit. 


14  . 


6.  DISPERSION  AT  HIGHER  FREC'w’ENCIES 

To  study  the  dispersion  at  resonant  and  igner  frequencies,  we 
must  consider  higher  powers  of  ka  ,  and  this  implies  that  a  large 
number  of  terms  !Y  ^  ,  Z  ^  )  musc  be  <epc  m  the  expansion  of  the 
average  field.  This  is  best  done  numerically.  A  block  diagram  of 
the  FORTRAN  program  written  for  this  purpose  is  shown  in  Fig.  1. 

The  blocks  identify  major  subroutines  which  perform  the  foliowma 
functions : 


MAIN  The  main  program  sets  up  the  basic  loops  so  as  to  calculate 
K  for  various  frequencies  (ka)  and  concentrations  ic) 

RDDATA  This  subroutine  is  used  to  input  data,  e.g.,  scatterer  sice 
;ka) ,  concentration  (c),  matrix  sices,  etc. 

RTINT  This  subroutine  calculates  the  initial  guess  for  X  in  the 
Rayleigh  limit  at  a  given  concentration  c  using  Eqs.  .  3Ci 
or  (31)  . 


TMAT  This  subroutine  calculates  the  T-matrix  for  given  scatterer 
chape  and  sice  (ka)  .  Current  maximum  size  is  40  •  40. 

CGRTQ  This  subroutine  searches  for  the  root  in  the  complex  plane 
igiven  an  initial  guess)  by  attempting  to  force  the  deter¬ 
minant  of  a  coefficient  matrix  (C)  »  0. 


AUX 

AB 


This  subroutine  sets  up  the  coefficient  matrix  according  to 
Eqs.  (14a, bl  .  Maximum  sice  is  40  •  40. 


This  subroutine  calculates  the  factors 
in  Appendix  B. 


oo ' 
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as  given 


TRIXJ  This  subroutine  calculates  the  Wigner  3-j  coefficients. 


CXMTX  This  subroutine  calculates  the  det  C  for  a  given  K  using 
standard  Gauss  elimination. 


The  computational  procedure  is  based  on  forming  the  coefficient 
matrix  C  according  to  Eqs.  (24a, b).  For  a  given  ka  ,  the  roots 
of  the  equation  det  C  »  0  are  searched  m  the  complex  K-plar.e 
using  an  iterative  root  searching  algorithm  which  employs  Muller's 
method.  Good  initial  guesses  for  K  were  provided  in  the  Rayleigh 
limit  by  Eqs.  (30)  or  (31),  and  these  could  be  used  systematically 


Attenuation  Constant  [2/3  =  2  im  (K/k)] 
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Figure  3.  Block  diagram  of  the  computer  program 


Figure  4.  Attenuation  vs  concentration  for  a  random  distri¬ 
bution  of  dielectric  spheres  at  ka  ■  11.3 
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to  obtain  convergence  at  increasingly  higher  values  or  ka.  Simiiar- 
I / ,  the  dependence  of  K  on  concentration  at  a  fixed  frequencv 
could  be  computed.  The  real  part  of  K  ,  ,  determines  the 

prase  velocity  while  the  imaginary  part,  K,  determines  the 
coherent  attenuation. 


7.  COMPUTATIONS  AND  COMPARISONS  WITH  EXPERIMENTS 

A  major  aim  of  the  computational  method  presented  here  is  to  pro¬ 
vide  a  means  of  studying  the  dispersion  characteristics  of  discrete 
random  media  at  higher  scatterer  concentrations.  Very  few  laboratory 
measurements  of  coherent  wave  attenuation  and  phase  shift  can  be 
found  in  the  open  literature.  However,  results  of  two  such  experi¬ 
ments  are  presented  here  [13,14]  and  compared  with  computations. 

The  first  experiment  refers  to  the  work  of  Hawley  et  ai.  [13] 
who  measured  the  coherent  wave  attenuation  and  phase  shift  through  a 
random  assembly  of  dielectric  spheres  (•:  =  1.0  34  )  blown  about  by 

turbulence  producing  fans.  The  measurements  were  performed  at  fixed 
microwave  frequency  corresponding  to  ka  «  11.8  and  at  various  con¬ 
centration  levels  up  to  maximum  packing.-  The  attenuation  measure¬ 
ments  are  shown  in  Fig.  4  together  with  computations  using  the 
present  theory  and  that  used  by  Hawley  et  al.  The  present  theory 
agrees  well  with  measurements  up  to  c  »  0.32  ,  beyond  which  the 
attenuation  is  predicted  to  decrease  and  goes  to  zero  at  c  *  0.54. 
The  dashed  curve  is  based  on  the  concept  of  the  "effective"  concen¬ 
tration  which  accounts  for  the  decreasing  available  volume  as  the 
number  density  increases.  The  "effective"  concentration  equals 
c'(l-c)  and  the  dashed  curve  in  Fig.  4  reflects  the  attenuation 
constant  as  a  function  of  this  altered  concentration.  The  computa¬ 
tions  now  agree  well  with  the  measurements.  The  conerer.t  phase 
shift,  -  ,  relative  to  free  space  is  plotted  as  a  function  of  con¬ 
centration  c  m  Fig.  5  where  the  experimentally  adjusted  values 
and  the  computations  are  in  good  agreement.  Attenuation  m  the 
Rayleigh  limit  (ka  »  0.05 )  for  the  same  scatterer  properties  is 
shown  in  Fig.  6  based  on  the  present  theory.  The  computations  fail 
for  o  -  0.125  since  Im(K)  becomes  negative.  This  feature  is 
repeated  in  the  Rayleigh  limit  for  other  values  of  as  shown 

later.  At  very  low  concentrations,  the  form  of  the  pair  correlation 
function  as  given  in  Eq.  (17)  seems  valid.  As  the  concentration 
increases,  the  available  volume  for  the  other  scatterers  decreases. 
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Figure  5.  Phase  shift  vs  concentration  or  number  of  spheres 
in  sample  volume  for  a  random  distribution  of  dielectric 
spheres  at  ka  *  11.8 


Figure  6.  Attenuation  factor  vs  concentration  for  a  random 
distribution  of  dielectric  spheres 
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Assuming  that  the  centers  of  04-1)  scatterers  are  randomly  located 
in  a  giver,  volume  V,  the  minimum  available  volume  for  locating 

the  center  of  the  N-th  scatterer  is  7^  =  v  -  i  -(2a)^!N-ll  *  V(l-3c' 

where  it  is  assumed  that  each  scatterer  center  is  surrounded  by  a 
hole  of  radius  2a  and  that  these  holes  do  not  interpenetrate.  Now 
V  -  0  as  c  -  0.125  so  this  may  explain  why  the  computations 
cause  lm(X)  •  0  for  c  •  0.125  ,  at  least  in  the  Rayleigh  limit. 
Figure  7  shows  computations  in  the  Rayleigh  limit  for  spheres  with 
=  3.158  and  ka  =  0.05  compared  with  the  analytical  results  of 
Twerskv  [2]  who  obtained  the  leading  effects  of  pair -correlation  at 

low  frequencies.  His  formula  reduces  to  Im(K)  »  r.  :  W  where 
4  2  2  o  s 

W  =  (1-c)  /(l+2c)  ,  n_  is  the  number  density  and  :  is  the 

o  J  s 

total  scattering  cross  section  of  a  lossless  sphere.  The  agreement 
is  good  for  c  <  0.05  while  great  discrepancy  is  exhibited  at 
higher  concentrations.  The  concept  of  the  “effective"  concentration 
does  not  significantly  alter  the  results  at  such  low  concentrations. 
To  obtain  better  agreement  at  higher  concentrations  it  is  necessary 
to  use  a  more  realistic  form  for  the  pair-correlation  function  de¬ 
fined  in  Ea.  (17).  A  factor  containing  some  dependence  of 
p(r^  r  )  on  r  .  has  to  be  introduced,  but  the  form  for  this 

i  •  *  J 

dependence  constitutes  a  difficult  problem  in  statistical  geometry 
(181  . 

A  second  set  of  coherent  wave  measurements  by  Olsen  and  Kharadiy 
(14]  at  relatively  low  concentrations  ic  -  0.007,  0.0141  and  at  a 
single  microwave  frequency  (ka  *  4.67)  was  available  for  comparison 
with  computations.  Their  measurement  procedure,  based  on  ensemble 
averaging  made  on  a  random  collection  of  dielectric  spheres 
(-:  *  2.26),  exhibited  greater  control  on  the  statistics  of  the 

scatterer  distribution  yielding  accurate  measurements  with  low  stan¬ 
dard  errors.  Attenuation  as  a  function  of  ka  is  shown  in  Fig.  3 
at  c  *  0.014.  The  experimental  value  is  also  shown  with  its  esti¬ 
mated  standard  error.  Figure  9  shows  the  same  results  at  c  *  3. OCT. 
As  expected,  the  agreement  between  theory  and  experiment  is  excellent 
at  these  low  concentrations.  Figures  10  and  11  show  the  phase  shift 
relative  to  free  space  as  a  function  of  ka  for  c  *  0.014  and 
o.on7,  respectively,  together  with  the  experimental  values.  Again 
agreement  is  good  and  within  the  estimated  standard  errors.  Figure 
12  shows  the  computation  of  attenuation  vs.  concentration  at 
ka  =  4.67.  The  values  based  on  single  scattering  theory,  see 
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figure  7.  Attenuation  factor  vs  concentration  for  a  random 
distribution  of  dielectric  spheres 


Figure  3.  Coherent  field  amplitude  vs  ka  for  a  random  dis¬ 
tribution  of  dielectric  spheres:  comparison  with  experimenta 
value  of  Olsen  and  Kharadly  (1975)  at  ka  »  4.67 
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Figure  II.  Phase  shirt  vs  ka  for  a  random  distribution  of 
dielectric  spheres;  comparison  with  experimental  value  of 
Olsen  and  Kharacly  (1975)  at  ka  =  4.67 
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Figure  12.  Attenuation  vs  concentration  for  i  random  dist 
bution  of  dielectric  spheres 
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iq.  i  3 4 ^  ,  are  also  shown  together  with  the  dasr.ed  curve  oased  an  the 
"effective"  concentration  concept.  In  this  case,  single  scattering 
theory  seems  to  be  -,'alid  for  c  3.35.  The  computations  cause 
I.T!  y :  •  3  for  c  _  0.6  wmle  the  dasr.ed  curve  does  not  exmcit 

this  phenomenon.  Again,  the  validity  of  the  compute t von s  at  higr.er 
concentrations  is  not  established .  Previous  computations  presentee 
m  ;6,7]  for  concentrations  of  0.20,  as  a  function  of  ka  or 
frequency) ,  exhibit  certain  "null1'  characteristics  fer  the  attenua¬ 
tion  at  certain  ka  values.  This  phenomenon  is  not  real  and  is 
caused  by  Im ; X )  0  at  these  values.  However,  one  computations 

presented  m  (6,7]  for  c  *  3.05  and  3.13  are  correct. 

3.  QUASI -CRYSTALLINE  APPROXIMATION  <CCA 

It  is  instructive  to  examine  the  physical  implications  of  terms 

of  the  form  <  b' _ -  that  occur  in  Ec .  -21).  This  is  tr.e  excitmc 

-  mp  --  3 

field  coefficient  of  the  3-th  scatterer  when  the  positions  of  all 
scatterers  except  the  3-th  and  any  other  scatterer  denoted  by  *  1' 
are  averaged  over.  When  the  numbei  of  particles  m  the  system  13 
large  this  may  be  thought  of  as  the  field  exciting  the  i-th  scat¬ 
terer  in  an  effective,  macroscopically  homogeneous  medium  containing 
two  sc3tterers  '1’  and  ' 3  * .  The  QCA  implicitly  omits  all  multiple 
scattering  processes  that  can  take  place  between  ' i'  and 
Efforts  have  been  made  to  restore  such  scattering  processes  by 
Twerskv  [2]  and  by  Schwartz  and  Ehrenreich  [19]  who,  in  a  different 
context,  discuss  the  contribution  of  clusters  of  two  or  more 
particles  as  in  systems  with  short  range  order. 

In  the  classical  context,  it  is  also  important  to  discuss  the 
dependence  of  QCA  on  the  frequencies  under  consideration  as  well  as 
on  the  concentration  of  the  scatterers.  Scattering  in  a  two  scat¬ 
terer  configuration  is  quite  dependent  on  the  ratio  of  the  distance 
’Id'  between  the  scatterers  to  the  overall  dimensions  '2a'  cf  tr.e 
scatterer.  Numerical  computations  suggest  f.nat  multiple  scattering 
effects  at  any  frequency  of  the  incident  wave  will  be  small  if  d,  a 
is  large.  Thus,  the  2CA  is  expected  to  be  good  for  sparse  concen¬ 
trations,  •'  1%.  At  such  concentre tions  the  pair  correlation  func¬ 
tion  as  used  in  Eq.  17  is  also  expected  to  be  good.  The  analytical 
resuits  we  obtain  at  long  wavelengths  give  ample  proof  of  this  as 
shown  m  the  previous  sections.  Most  authors  cited  previously  as 
well  as  Talbot  and  Willis  (20]  in  their  recent  work  seem  to  agree 
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that  the  QCA  is  good  at  sparse  concentrations .  However,  most  of 
these  papers  present  only  long  wavelengt.n  results.  Further,  lax  (21 
comments  lotted  lately  fallowing  the  definition  of  the  QCA  that  it  is 
good  for  dense  systems  and  exactly  valid  for  systems  with  a  crysta.- 
ime  structure.  These  are  contradictory  observations  about  the  same 
approximation  and  needs  to  be  studied  further. 

Our  computations  for  various  concentrations ,  but  more  important¬ 
ly  fcr  various  values  of  Xa  < the  non— dimensional  wavenumber)  sug¬ 
gest  that  even  at  low  values  of  Xa  3.35  ,  the  model  fails  if  the 
volume  concentration  c  exceeds  12.5%,  whereas  at  xa  values  2.1 
or  more,  the  CCA  leads  to  reasonable  results  for  t.ne  bulk  propaga¬ 
tion  constant  at  all  values  of  the  concentration  m  spite  of  the 
poor  model  used  for  the  pair  correlation  function.  Vie  expect  all 
types  cf  multiple  scattering  effects  including  cluster  effects  to  be 
important  at  concentrations  '-0.1  or  more,  but  it  should  be  noted 
that  the  QCA  type  approximation  neglects  only  repeated  scattering 
between  pairs  or  within  a  croup  of  scatterers.  It  is  also  important 
to  note  that  at  wavelengths  comparable  tc  obstacle  size  and  higher, 
i.e.,  Xa  •  3.0  ,  the  scattering  is  mostly  m  t.ne  forward  direction 
Thus,  m  this  case  repeated  scattering  should  not  be  important, 
since  the  backscattered  wave  is  sigmf icaneiy  smaller  than  the  for¬ 
ward  scattered  wave.  This  would  help  satisfy  the  QCA  and  may  ex¬ 
plain  why  the  computations  shown  m  Fig.  4  for  Xa  »  11.3  are  m 
reasonable  agreement  with  experiment  at  all  values  of  the  effective 
concentration.  It  would  appear  that  in  the  context  of  classical 
systems,  Lax’s  statement  about  the  validity  of  the  QCA  for  dense 
systems  should  be  qualified  by  the  phrase,  ’at  high  frequencies.’ 

It  would  be  interesting  to  study  how  tne  results  will  change  with 
improved  models  of  the  pair  correlation  function.  Then  it  wou»d  be 
possible  to  comment  on  the  sensitivity  of  the  results  to  the  effect 
of  the  pair  correlation  function  and  the  QCA  separately . 

9.  RECOMMENDATIONS  FOR  FUTURE  WORK 

It  is  obvious  from  the  preceding  discussions  on  the  QCA  as  well 
as  the  numerical  results  that  the  two  major  improvements  required 
are  for  the  QCA  as  well  as  the  pair  correlation  function,  so  that 
good  results  can  be  obtained  for  all  concentrations  even  at  long  and 
intermediate  values  of  the  wavelength.  In  a  review  article.  Lax 
(21)  has  suggested  that  in  the  quantum  mechanical  context,  the  QCA 
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could  oe  improved  by  using  modified  propagators  for  t he  fields.  In 
the  classical  context,  this  implies  that  on  the  average,  single  par¬ 
ticle  scattering  taxes  place  m  a  macroscopicaliy  homogeneous  medi¬ 
um,  and,  in  this  respect,  this  idea  is  the  sane  as  the  cone rent 
potential  approximation  ,CPA;  of  Solid  State  Physics.  The  repeated 
multiple  scattering  between  pairs  of  scatterers  or  cluster  effects 
can  be  improved  by  maxing  the  self  consistent  approximation  -.SCA)  m 
addition  to  the  CPA. 

For  tne  purpose  of  discussion  of  these  ideas  within  the  T -matrix 
formalism  civen  earlier,  we  denote  by  E^  and  E,  the  r le Ids  scat- 
tered  by  and  exciting  the  j-th  scat terer,  respectively.  The  expan¬ 
sion  coefficients  of  these  fields  as  given  in  Eqs  .  (5,6>  are  denoted 

by  3 3  and  b-  ,  respectively,  omitting  all  subscripts. 

The  CPA  car.  be  expressed  succinctly  as 

33',  =  TC<!  b3^  -15! 

where  the  T-natrix  relating  the  exciting  and  scattered  field  coeffi¬ 
cients  is  evaluated  using  the  bulk  propagation  constant  X  for  the 
embedding  medium.  Thus  the  CPA  implies  that  the  field  scattered  by 
a  single  obstacle  ir.  the  presence  of  several  others  when  averaged 
over  the  position  of  all  scatterers  is  the  same  as  the  field  that 
would  be  produced  by  a  single  particle  embedded  in  a  microscopically 
homogeneous  medium  described  by  the  propagation  constant  K.  The  in¬ 
corporation  of  the  CPA  into  the  previous  formalism  involves  changes 
only  in  the  computations  and  a  redefinition  of  the  T-matnx.  It 
would  be  interesting  to  see  the  change,  if  any,  in  tne  numerical 
computations  as  a  result  of  invoking  the  CPA. 

The  idea  beni.nd  the  'self  consistent  approximation'  -SCA)  is 
somewhat  more  subtle.  From  the  discussion  in  the  section  on  the  2CA, 
it  is  now  clear  that  QCA-CPA  neglects  multiple  scattering  between 
two  fixed  scatterers.  The  SCA  as  defined  by  Schwartz  and  Ehrenreich 
(19!  restores  this  by  stating  that 

'B3> .  .  -  T (K)  <b3>  .  v  36 ) 
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where  f(X)  is  the  'T-matnx'  of  3catterer  '  j  ’  in  the  presence  of 
scatterer  'i',  in  the  effective  medium  with  propagation  constant  X. 
Expressions  for  T(K)  as  given  by  Peterson  and  Strom  [16]  may  be 
written  as 
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r  -  X)  »  3!?..  2  Til  -  ;  '  -r.  .  ,T 

:i  *  :  -r. . iTR ;r. . ! )  3i-r  2)  ■  ■* 

where  -  is  a  compact  notation  tor  the  translation  matrices  A  and 
3  introduced  in  Eq.  10).  The  R  matrix  is  simply  the  part  of  : 
that  is  regular  at  the  origin,  i.e.,  for  r  3  0.  All  matrices 
m  tc.  :3~>  are  obtained  using  the  buix  propagation  constant  for  tr.e 
host  medium. 

We  observe  that  T  K)  explicitly  depends  on  r, .  the  distance 

A  i 

between  ’  i’  ar.d  ’  ] 1  and  hence  will  be  involved  in  the  integrations 
m  Eqs.  A1-A12)  .  The  integration  procedure  will  r.c  longer  be  sim¬ 
ple  as  before  and  the  SCA  may  oe  rather  difficult  to  enforce  m  com¬ 
putations,  especially  if  more  realistic  models  are  chosen  for  tr.e 
pair  correlation  function. 

Improvements  to  the  pair  correlation  function  must  taice  info 
account  the  increase  m  short  range  order  with  increasing  concentra¬ 
tion  in  addition  to  enforcing  tr.e  condition  of  no  interpenetration 
of  particles.  Talbot  and  Willis  [20]  i r.  their  recent  wort  nave  re¬ 
viewed  several  models  that  depend  on  concentration.  They  discuss 
in  particular  two  models  by  Matern  .26'.  that  depend  on  the  concept 
of  an  available  volume  and  hence  are  valid  only  for  c  12.  :i. 
Talbot  and  Willis  also  suggest  the  Percus-VevicX  ,'2'S  ncdei  for  the 
pair  correlation  function  wmch  is  valid  for  o  •  0.3. 

Incorporation  of  the  CPA  as  well  as  improved  models  of  the  pair 
correlation  function  into  our  computations  are  m  progress.  We  nope 
that  they  will  shed  some  light  or,  the  sensitivity  of  multiple  scat¬ 
tering  theories  to  approximations  like  CCA  ar.d  -CA  as  a  function  of 
frecuency  and  scatterer  concentration.  Needless  to  say  additional 
experimental  resuits  are  required  for  comparison  with  these  compu¬ 
tations  . 
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Consider  the  following  integral  whic.n  appears  in  Ea.  23a): 


Y.1P  Aoin’  3 
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Use  of  Green's  cheorem  in  (AS)  reduces  it  to 
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where  Sx  refers  to  the  surface  of  a  sphere  of  larce  radius  cen¬ 
tered  around  the  j-th  soatterer  and  S-,^  refers  to  the  surface  of 

the  hole  centered  around  the  3-th  scatterer.  The  integral  over  S, 

c  a 

can  be  evaluated  m  closed  form 
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where  n  is  the  unit  outward  normal  to  S,a-  The  above  integral 
reduces  to 
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Upon  substituting  the  previously  defined  expressions  for  i  a.-.d 
3n  and  using  the  orthogonality  of  the  Legendre  polynomials  in  A' 
we  get 

-1—-X  (JH)*  ,y  „  ;  (JH)*  *  '  JH)  a  A3 

k2-X2  '  '2'*1’  'n 

where  (JH)  is  given  by  Eq  .  (25)  and  •  is  the  Kornecker  delta. 

Hence,  I  defined  in  A6)  reduces  to 
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wnere  I  _  refers  to  me  integration  over  3  wr.icr.  is  octai-ted  cy 
jsmcr  me  asvnptotic  expansion  for  h  kr^,  35  r  .  *  •'  on  5 j 

?!-.*  integral  germed  m  A 2'  now  reduce-;  to 
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Final  lv.  me  m  tear  a  *  given  m  A-  reduces  to 
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In  a  similar  manner ,  t he  renaming  integrals  . n  Eg.  11a. o  ran  oe 
evaluated . 
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The  trarslat  icna.  operators  .  n,n‘.  and  r. ,  r.  *  ,  are 
g i ven  l n  .  1  e  .  : 
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'be  above  f 

actors  can  be 

expressed 

in  terms 

of  factors  given  by 

'ruzan  [17] 

as  follows 
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ABSTRACT 


The  coherent  electromagnetic  wave  attenuation  in  an  infinite  medium 
composed  of  a  random  distribution  of  identical,  finite  scatterers  is 
studied.  A  self-consistent  multiple  scattering  theory  using  the  T-matrix 
of  a  single  scatterer  and  a  suitable  averaging  technique  is  employed.  The 
statistical  nature  of  the  position  of  scatterers  is  accounted  for  by  ensemble 
averaging.  This  results  in  a  hierarchy  of  equations  relating  the  different 
orders  of  correlations  between  the  scatterers.  Lax's  quasicrystalline 
approximation  (QCA)  is  used  to  truncate  the  hierarchy  enabling  passage  to 
a  homogeneous  continuum  whose  bulk  propagation  characteristics  such  as 
phase  velocity  and  coherent  wave  attenuation  can  then  be  studied.  Three 
models  for  the  pair  correlation  function  are  considered.  The  Matern 
model  and  the  well  stirred  approximation  (WSA)  are  good  only  for  sparse 
concentrations,  while  the  Percus-Yevick  approximation  (P-YA)  is  good  for 
a  wider  range  of  concentration.  The  results  obtained  using  these  models 
are  compared  with  the  available  experimental  results  for  dielectric 
scatterers  embedded  in  another  dielectric  medium.  Practical  applications 
of  this  study  include  radar  meterologv  and  communications  through 
hydrometers,  dust,  vegetation,  etc. 
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1 .  INTRODUCTION 


We  consider  the  propagation  of  plane  coherent  electromagnetic  waves 
in  an  infinite  medium  containing  identical,  lossless  randomly  distributed 
particles.  Our  aim  is  to  characterize  the  random  medium  by  an  effective 
complex  wave  number  K  which  would  be  a  function  of  the  particle  concentration, 
electrical  size  and  the  statistical  description  of  the  random  positions  of 
the  scatterers.  The  imaginary  part  of  K  describes  the  coherent  attenuation 
which  is  due  to  multiple  scattering  only  since  the  particles  themselves 
are  assumed  to  be  lossless.  The  understanding  of  the  behavior  of  Im(K)  as 
a  function  of  particle  concentration  c  and/or  frequency  ka  is  very 
important  in  many  practical  applications,  including  wave  propagation  in 
the  atmosphere  and  oceans  and  whenever  distribution  of  random  scatterers 
influence  electromagnetic  wave  behavior. 

The  theoretical  formulation  presented  here  closely  follows  the 
procedure  described  in  Varadan  et.  al.  [1979]  and  Bringi  et.  al.  [1981]. 

This  approach  is  based  on  a  self-consistent  multiple  scattering  theory 
and  relies  on  the  T-matrix  [Waterman  1971]  which  relates  the  field 
scattered  by  a  particle  to  an  arbitrary  exciting  field.  The  statistical 
description  of  the  random  position  of  the  scatterers  is  used  to  define  a 
configurational  average  which  results  in  a  hierarchy  of  equations  relating 
the  different  orders  of  correlations  between  the  scatterers.  Lax's 
[1952]  quas i -crys ta 1 1 ine  approximation  is  used  to  truncate  the  hierarchy 
which  results  in  the  usual  "hole-correction"  integrals.  Following  Twersky 
[1977,  1978  a,b],  a  radially  symmetric  pair-correlation  function  is 
introduced  and  approximate  models  are  chosen  from  Talbot  and  Willis  [1980]. 


The  "we  1 1 -st irred"  approximation  (WSA)  was  used  previously  by  Varadan 
et .  al.  [1979]  and  Bringi  et.  al.  [1981]  which  assumes  no  correlation 
between  the  particles  except  that  they  should  not  inter-penetrate .  In 
particular,  the  WSA  gives  unphysical  results  for  c  >_  0.125  at  the  Raleigh 
or  low  frequency  limit. 

In  this  paper,  we  consider  two  other  pair-correlation  functions, 
viz.  (i)  the  Matern  [1960]  model  and  (ii)  the  Percus-Yevick  [1957]  model 
for  a  classical  system  of  hard  spheres.  Computations  of  Im(K)  are  presented 
for  dielectric  scatterers  in  a  dielectric  medium,  using  the  above  three 
models  as  a  function  of  frequency  and  concentration.  We  also  compare  our 
solution  to  some  recent  optical  propagation  experiments  conducted  by 
Is'nimaru  [1981].  Sample  computations  are  also  presented  comparing  the  WSA 
and  the  single  scattering  approximation  for  a  rain  medium. 

2.  FORMULATION  OF  THE  PROBLEM 

Consider  N  identical,  finite  dielectric  scatterers  that  are  randomly- 
distributed  either  in  free  space  or  in  a  different  dielectric  medium.  The 
scatterers  are  homogeneous  with  a  relative  dielectric  constant  of  c^,  their 
centers  being  denoted  by  0^,  0,.  CL,  ...,  0^.  They  are  assumed  to  be  bodies 
of  revolution  with  symmetry  axis  parallel  to  the  z-direction.  Monochromatic 
plane  coherent  electromagntic  wave  is  assumed  to  propagate  along  the  symmetry 
axis  of  the  scatters  to  satisfy  the  condition  that  the  effective  medium  be 
isotropic  and  polarization  insensitive.  The  time  dependence  of  the  incident 
field  and  hence  the  fields  scattered  by  the  individual  scatterers  is  all  of 
the  form  exp (-jut)  and  this  is  suppressed  in  the  equations  that  follow. 
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Even  though  the  theory  presented  here  is  valid  for  spheroidal  scatterers 
(Varadan,  et.  al .  1981],  we  present  numerical  results  only  for  spherical 
scatterers  in  order  to  compare  our  results  with  available  experiments. 

Let  E^(r)  be  the  electric  field  arising  from  the  incident  plane  wave 
and  ET’(r)  the  field  scattered  by  the  i-th  scatterer.  Both  these  fields 
satisfy  the  vector  Helmholtz  equation.  The  problem  at  hand  reduces  to 
computing  the  total  wave  field  at  any  point  outside  the  scatterers, 
satisfying  the  appropriate  boundary  condition  on  the  surface  of  the  scatterer 
and  radiation  conditions  at  infinity. 

The  total  field  at  any  point  outside  the  scatterers  can  be  interpreted 
as  the  sum  of  the  incident  field  and  the  fields  scattered  by  all  the 
scatterers,  which  can  be  written  as 


where  Eh(p^)  is  the  field  scattered  by  the  i-th  scatterer  at  the  observation 
point  r.  However,  the  field  that  excites  the  i-th  scatterer  is  the  incident 
field  E°  plus  the  fields  scattered  from  all  other  scatterers  except  the  i-th. 
The  term  exciting  field  E  is  used  to  distinguish  between  the  field  actually 
incident  on  a  scatterer  and  the  external  incident  field  E°  produced  by  a 
source  at  infinity.  Thus,  at  a  point  r  in  the  vicinity  of  the  i-th 
scatterer,  we  write 


a  < 


<  2a 
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where  'a'  is  a  typical  dimension  of  the  scatterer. 

The  exciting  and  scattered  fields  for  each  scatterer  can  be  expanded 
in  terms  of  vector  spherical  functions  with  respect  to  an  origin  at  the 
center  of  that  scatterer: 
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E*(r) 


f  F  F  b1.  Re  y  ,  (?.) 

Tin a  t tnc1  U 

c= 1  i=l  n=0  o=l 


(3) 


=  7  b  Re  ^ 

in  rn 
in 


-+-S  —  i 

E  (r)  =  )  B 

l  “  rn 

rn 


where  the  vector  spherical 
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functions  are  defined  as 

hF)»r>i 
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In  equations  (3-6),  k  is  the  wave  number;  h£*^  is  the  Hankel  function  of  the 
first  kind  and  the  Y,  (6,4)  are  the  normalized  spherical  harmonics  defined 
with  real  angular  functions.  In  Equation(3),  the  exciting  field  is  expanded 
in  terms  of  the  regular  (Re)  basis  set  (Re  ^  )  obtained  by  replacing  h^^  in 
Equations  (5-6)  by  j  ,  the  spherical  Bessel  functions  of  the  first  kind.  Thus, 
the  choice  of  the  basis  set  in  Equation  (4)  satisfies  the  radiation  condition 
at  infinity  for  the  scattered  field,  while  the  choice  in  (3)  satisfies  the 
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regular  behavior  of  the  exciting  field  in  the  region  a  <  |p^|  <  2a.  The 
superscript  i  on  the  basis  functions  refer  to  expansions  with  respect  to 
0.,  and  bL  and  81  are  the  unknown  exciting  and  scattered  field  coefficients. 
We  also  expand  the  incident  field  in  terms  of  vector  spherical  functions: 


I  (r)  =  e 


ikz  •  r 
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tn 


a  Re  'p  (e . ) 
tn  in  i 
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where  a  are  the  known  incident  field  coefficients, 
tn 

The  unknown  coefficients  b1  can  be  related  to  B1  by  means  of  anv 

tn  in 

convenient  scattering  operator,  in  this  case  we  employ  the  T-matrix  as  defined 
by  Waterman  [1971]: 


B1  =  l  T1  ,  ,  b1,  , 
tn  V  ,  rn.T'n'  t ’n' 
t'n' 


Substituting  Equations  (5),  (4)  and  (7)  in  (2),  we  obtain 
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l  b1  Re  l1  =  elkz’ri  l  a  Re  ^  \  l  BJ  r1 
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Since  the  field  quantities  are  expanded  with  respect  to  centers  of  each 
scatterer,  we  obtain  Equation  (9)  with  basis  functions  expanded  with  respect 
to  i-th  and  j-th  centers.  In  order  to  express  them  with  respect  to  a  common 
origin  Ch  ,  we  employ  the  translation  and  addition  theorems  for  the  vector 
spherical  functions  [see,  for  example,  Bostrom,  19S0] 

which  may  be  written  in  a  compact  form  as  follows 


= 
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where  P..  =  r.-r.  is  the  vector  connecting  0.  to  0. ,  o  ,  , 
ij  i  j  j  i  xn.x'n' 


is  the 


translation  matrix  for  the  vector  functions  and  R  .  , 

tn ,  t  n 


is  a  matrix  with 


spherical  Hankel  functions  in  aTnitini  replaced  by  spherical  Bessel  functions. 

Employing  Equations  (8)  and  (10)  in  (9)  and  using  the  orthogonality 
of  the  vector  spherical  basis  functions,  we  obtain  the  following  set  of 
coupled  algebraic  equations  for  the  exciting  field  coefficients  b*  : 

1  =  eik:  •  rj_  +  y  y  y  a  (p.  ,)Tj  ,  ,  bj ,  ,  (11) 

m  tn  .4,  L  r  ,  xn.x’n'^  ij  Tn.i’n'  x'n' 
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From  Equation  (11),  it  can  be  seen  that  the  exciting  field  coefficients 

of  the  i-th  scatterer  explicitly  depend  on  the  position  and  orientation  of 

the  other  scatterers.  In  this  paper,  we  consider  a  random  distribution  of 

spherical  scatterers  and  the  case  when  N  -*■  00  and  the  volume  occupied  by  the 

scatterers  V  -*■  »  such  that  N/V  =  n  is  a  finite  number  densitv.  For  such 

o 

distribution,  a  configurational  average  of  Equation  (11)  can  be  made  over 

the  positions  of  all  scatterers  [see  Varadan  et .  al.,  1981]  with  QCA 

[Lax,  1952]  to  arrive  at  an  equation  for  the  configurational  average 

<b1  >.  of  the  exciting  field  coefficients  with  one  scatterer  fixed: 
xn  i 
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where  p(r^|r^)  is  the  two  particle  joint 
the  above  equation,  we  have  assumed  that 
We  now  assume  that  the  average  field 
propagates  in  a  medium  with  an  effective 
in  the  direction  of  the  original  incident 
medium : 


probability  density.  In  obtainin 

all  the  scatterers  are  identical. 

<bX  >.  (the  coherent  field) 
tn  i 

complex  wave  number  K  =  (K^+iK-,)z 
field  in  the  discrete  random 
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Substituting  Equations  (13)  and  (14)  in  Equation  (12)  and  invoking  the 
extinction  theorem  to  cancel  the  incident  wave  term  in  (12),  we  obtain 
the  following  eqqations  for  the  unknown  amplitudes  and 
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In  the  aoO'/e  equation,  c  =  4n  n^a  /3  is  the  effective  spherical  concentration. 

For  plane  waves  propagating  parallel  to  the  rotational  axis  of  symmetry  of 
scatterers,  only  1=1  contributes,  and  also  only  certain  combinations  of  n 
yield  non-zero  1 -matrix  elements  which  are  used  in  Equations  (IS)  and  (16). 

In  Equation  (17),  g(.x)  is  the  pair  correlation  function  which  depends  only 
on  jxj  =  jp\j  |  due  to  translational  invariance  of  the  system  under  consideration. 
To  obtain  expressions  for  g(x),  a  description  of  the  interparticle  forces  is 
needed.  In  our  statistics,  the  dielectric  scatterers  are  assumed  to  behave 
like  effective  hard  spheres  of  radius  'a'  where  'a'  is  the  radius  of  the 
circumscribing  sphere,  see  Figure  1.  Wertheim  [1965]  has  obtained  a  series 
solution  of  the  integral  equation  for  the  pair  correlation  function  derived  by 
Percus  and  Yevick  [1958]  for  an  ensemble  of  hard  spheres.  Throop  and  Bearman 
[1965]  have  used  the  Wertheim  result  and  provided  tabulated  values  of  g(.\) 
as  a  function  of  x  for  several  values  of  c.  Plots  of  g(x)  vs  x  is  shown  in 
Figure  2. 

At  low  values  of  concentration  c,  g(x)  ~  1,  see  Figure  2  and  hence  the 
integral  in  Equation  (17)  is  negligible  which  results  in  a  system  of  uncorrelated 
hard  particles.  This  is  what  has  been  referred  to  as  the  well  stirred 
approximation  (tvSA)  and  yields  the  'hole  correction  integral'  as  outlined  by 
Fikioris  Waterman  [196J]  and  by  us  earlier.  If  g(x)  >  1,  one  can  regard  the 


10 


Equation  (17)  as  a  modified  'hole  correction  integral'  which  is  of  the  same 

form  as  used  by  Twersky  [1977,  1978]. 

Equations  (.15)  and  (16)  are  simultaneous  linear  homogeneous  equations  for 

the  unknown  amplitudes  Y  ..  For  a  nontrivial  solution,  we  require  that 
r  tanu 

the  determinant  of  the  truncated  coefficient  matrix  C  vanishes,  which  yields 
an  equation  for  the  effective  wave  number  K  =  (K^  +  iK.,)  in  terms  of  k  and  the 
T-matrix  of  a  scatterer.  This  is  the  dispersion  relation  for  the  scatterer 
filled  medium.  The  real  part  of  K  relates  to  the  phase  velocity  while  the 
imaginary  part  relates  to  coherent  attenuation  in  the  medium. 

3.  NUMERICAL  COMPUTATIONS 

In  the  low  concentration  limit,  c  ->•  0,  it  is  well  known  that  the  single 
scattering  approximation  (SSA)  is  valid  so  that  Im(K/k)  is  given  by 

3  ^ext 

Im(K/k)  =  J  c  -gi-  (20) 

-> 

where  Qext  is  the  normalized  (with  respect  to  Tia“)  extinction  cross  section 

of  a  sphere  of  radius  'a'.  An  important  problem  is  propagation  in  a  rain 

medium  where  the  single  scattering  approximation  has  been  widely  used.  Indeed, 

even  under  very  heavy  rain,  the  concentration  rarely  exceeds  0.01  and  is 
-4 

typically  around  10  .  We  have  compared  our  theory  using  WSA  with  Equation 

(20)  for  a  distribution  of  spherical  water  drops  of  radius  0.1  cm  with  ka  in 
the  range  0.1  <_  ka  <_  3.  The  refractive  index,  which  is  a  function  of  frequency, 
is  taken  from  Ray  [1972].  In  Figure  (3),  we  show  the  attenuation  constant  y 

-3 

defined  as  4x  Im(K)/Re(K)  as  a  function  of  ka  using  the  WSA  for  c  =  10  ”,  10  , 

-4 

and  10  which  is  to  be  compared  with  Figure  (4)  which  uses  SSA.  We  note  that 


both  solutions  yield  nearly  identical  results.  In  Figure  (5),  we  show 
computation  of  y  vs  concentration  for  different  ka  values  using  the  WSA. 
Again  the  SSA  is  seemed  to  be  excellent  for  the  rain  medium. 


We  now  present  computations  for  a  random  medium  model  used  by  Ishimaru 
[1981]  for  the  optical  propagation  experiments.  The  scatterers  are  latex 
spheres  of  diameter  0.107y  immersed  in  water  with  incident  wavelength 
X  =  0.6u.  In  the  Rayleigh  limit,  Twersky  [1978b]  has  given  an  expression  for 
Im(K/k)  by  considering  the  leading  effects  of  the  pair-correlation: 

2 

W  (21) 

where  is  the  relative  dielectric  constant  and  IV  is  the  packing  factor 
given  by 


Im(K/k)  =  c(ka)‘ 


£  -1 
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W  =  l1'-)  -=■  =  1  +  24c  /  x2[g(x)-l]  dx  .  (22) 

(l+2c)  o 

In  Figure  (6),  we  show  Im(K/k)  as  a  function  of  concentration  c  using  Equation 
and  the  present  theory  employing  the  WSA,  the  P-YA  and  the  Matern  model.  The 
Matern  [1960]  model  is  completely  analytic  and  is  valid  for  c  <  0.125.  We  note 
that  Equation  (22)  and  the  P-YA  are  identical  while  both  the  Matern  model  and 
the  WSA  fail  for  csA0.04,  and  in  fact  they  give  unphysical  results  for 
c  >  0.12S. 

In  Figure  (7),  we  show  the  comparison  between  the  computation  and  the 
two  measured  values  at  c  =  0.01  and  0.10  given  by  Ishimaru  [  1 9 S 1  ]  .  We  note 
that  the  ka  value  is  0.56  and  that  multiple  scattering  effects  are  seemed  to 
be  important  even  at  c  =  0.01.  The  measured  values  at  c  =  0.01  and  0.1  are  in 
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very  good  agreement  with  both  the  WSA  and  P-YA  while  the  SSA  consistently 
overestimates  the  effective  coherent  attenuation.  Also,  for  c  >  0.10  where 
measurements  are  not  available  at  the  present  time,  we  feel  that  only  the 
P-YA  predicts  the  correct  behavior  of  Im(K/k).  In  Figure  (8),  we  show  the 
variation  of  Im(K/k)  with  ka  for  c  =  0.21  and  compare  the  results  using  the 
SSA,  the  WSA  and  the  P-YA.  Values  for  the  WSA  for  ka  <_  0.7S  are  not  shown 
since  the  solution  fails  [Im(K/k)  <  0]  in  this  region.  However,  as  ka 
increases  it  appears  that  the  WSA  tends  to  merge  with  P-YA  for  ka  ^  3.0. 

The  SSA  on  the  other  hand  predicts  a  higher  attenuation  than  either  the  WSA 
or  the  P-Ya. 
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Abstract 


A  scattering  matrix  theory  is  presented  tor  studying  the  multiple 
scattering  of  both  longitudinal  and  transverse  elastic  waves  in  a  medium 
containing  a  random  distribution  of  inclusions  or  voids  of  arbitrary  shaoe. 

A  statistical  analysis  with  QCA  and  Percus - Yevick  pair  correlation  function 
isthen  employed  to  obtain  expressions  for  the  average  amplitudes  of  the 
coherent  fields  which  may  be  solved  to  yield  the  bulk  or  effective  properties 
of  the  inhomogeneous  medium.  Suggestions  for  incorporating  CPA  in  conjunction 
with  QCA  so  that  materials  with  dense  concentration  of  inclusions  can  be 
considered  are  also  given. 


Introduction 

In  recent  years,  considerable  effort  has  been  devoted  to  promoting 
the  development  of  elastomeric  absorbing  materials,  containing  a 
distribution  of  cavities  and  inclusions,  which  are  bonded  to  submerged 
structures  to  control  the  sound  radiated  by  these  structures  as  well  as 
to  modify  their  acoustic  reflection  characteristics  (echo  reductionl.  To 
use  such  absorbing  layers,  it  is  important  to  determine  how  their  physical 
properties  such  as  density,  thickness  and  effective  elastic  moduli,  and 
material  composition  such  as  distribution  and  orientation  of  the  inclusion 
and  their  sice  distributions  affect  the  acoustical  behavior  of  any  actual 
sturcture  coated  with  that  material. 

The  waves  incident  on  such  inhomogeneous  media  undergo  multiple- 
scattering  due  to  the  presence  of  inclusions  thus  reducing  the  scattering 
amplitude  or  cross  section  by  absorption  and  attenuation  of  waves.  The 
attenuation  depends  critically  on  the  material  properties  of  the  host 
medium  (matrix)  and  inclusions,  the  distribution  of  the  inclusions  and  the 
frequency  of  the  incident  wave.  The  problem  is  very  difficult  and  to  our 
knowledge,  rigorous  theories  with  numerical  results  are  not  available  in 
the  literature. 

In  multiple  scattering  theories,  approximations  are  usually  made  at 
a  very  early  stage  for  a)  the  geometry  of  the  inclusion,  b)  the  site  of 
the  inclusion  relative  to  the  wavelength  of  incident  wave,  ana  c) 
distribution  of  the  inclusions  in  the  matrix  median;.  The  approx  mat  ions 
with  respect  to  geometry  and  site  are  related.  If  t ho  inclusion  is  small 
compared  to  the  incident  wavelength,  it  is  not  possible  to  "see"  exact 
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details  of  the  inclusion  and  usually  one  is  content  to  obtain  the  gross 
scattering  properties  of  the  inhomogeneous  medium.  This  is  the  so-called 
Rayleigh  or  low  frequency  limit,  and  yields  corrections  to  the  solution 
for  point  scatterers.  As  far  as  the  distribution  of  the  inclusions  is 
concerned,  one  either  has  regular  arrays  of  inclusions  or  a  random 
distribution.  In  the  former  case,  one  performs  a  lattice  sum  while  in 
the  latter  case,  one  employs  a  configurational  averaging  procedure.  If 
the  concentration  of  inclusions  is  small,  i.e.,  the  inclusions  are  sparsely 
distributed,  we  may  use  a  single  scattering  or  first  Born  approximation. 
Approximations  have  been  employed  by  many  authors  and  the  correspond-n  ; 
effective  properties  of  the  medium  were  studied  at  the  low  frequencies  and 
low  concentrations,  see  for  example,  Waterman  and  Truell  [1],  Merkulova  [2], 
Chaban  [3,4],  Chatterjee  and  Mai  [S] ,  Domanv,  Gubernatis  and  Krumhansl  [6], 
Korringa  [7],  Kroner  [8],  Datta  [9]  and  the  references  therein.  Actually 
the  real  problem  warrants  a  rigorous  multiple  scattering  theory  and  a 
computational  approach  to  study  the  frequency  dependent  properties  of  the 
inhomogeneous  media  which  will  be  valid  for  frequencies  comparable  to 
scatterer  sice  and  for  a  wide  range  of  concentrations,  shapes  and  sices. 

Recently,  the  present  investigators  have  developed  a  multiple  scattering 
formalism  by  introducing  the  concept  of  a  T-matrix  for  individual  inclusions 
that  makes  the  formulation  more  general  and  applicable  to  a  variety  of 
different  scatterers,  see  Refs.  [10-20],  The  method  also  lends  inself  to 
numerical  computations  for  higher  frequencies  of  the  incident  plane  wave  as 
well  as  more  realistic  geometries  for  the  inhomogeneities.  The  dynamic 
elastic  properties  of  composite  elastic  media  have  been  studied  in  [20] 
using  tins  formulation,  and  the  concept  of  an  average  frequency  dependent 
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elastic  stiffness  tensor  following  the  work  of  Bedeaux  and  Mazur  [21], 
and  Varadan  and  Vezetti  [22].  The  results  seem  to  be  promising  for  future 
research  in  this  area.  In  Ref.  [20],  we  have  shown  that  a  Clausius -Mosott  i 
type  formula  for  the  average  shear  modulus  can  be  recovered  in  the  low 
frequency  limit.  For  higher  frequencies,  we  have  obtained  the  dynamic 
properties  for  a  range  of  frequencies.  The  extension  of  the  theory 
presented  in  [20]  to  acoustic  and  elastic  wave  scattering  will  be  useful 
for  Maval  applications. 

The  present  state  of  the  art  is  as  follows:  the  statistical  considerations 
seem  to  be  the  most  difficult  for  three  dimensional  inclusions  and  the  least 
amount  of  progress  has  been  made  in  this  area.  All  formalisms  that  involve 
ensemble  averaging  result  in  a  heirarchy  of  equations  for  the  average  fields 
that  involve  higher  and  higher  order  correlation  functions.  This  heirarchy 
must  be  truncated  in  some  fashion.  Foldy  [23]  approximated  the  field  incident 
on  a  scatterer  by  the  average  field  itself.  Lax  [24]  was  the  first  to  use 
a  quasi-crystalline  approximation  which  involves  the  two  particle  correlation 
function.  At  the  moment,  only  the  'hole  correction'  has  been  taken  into 
account  in  a  systematic  way.  Bose  and  Mai  [25]  have  tried  correlation 
functions  that  fall  off  exponentially  with  distance.  Recently,  Twersky  [2o] 
has  used  the  scaled  particle  equation  cf  state  of  a  gas  of  hard  spheres  to 
obtain  improvements  to  the  hole  correction  integral.  The  T -matrix  formalism 
employs  Lax's  quasicrystalline  approximation  (.QCA) ,  the  hole  correction 
integral  and  results  in  a  set  of  equations  that  must  be  solved  in  a  self- 
consistant  manner. 

In  this  paper,  a  radially  symmetric  pair-correlation  function  given  h> 
Percus-Yevick  (P-YA)  integral  equation  [2“]  is  introduced  which  gives 
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improvements  to  the  hole  correction  integral.  The  "we  1 1 -stirred"  approxi.mat ion 
(WSA)  was  used  previously  by  us  which  assumes  no  correlation  between  the 
scatterers  except  that  they  should  not  interpenetrate.  The  WSA  seems  to 
depend  on  concentrat ion  and  frequency.  At  low  frequency  or  Rayleigh  limit, 

WSA  gives  good  results  up  to  concentration,  c  <_  0.04  and  unphysical  results 
for  c  2.  0.12S  [28],  However,  at  higher  frequencies  and  higher  concentrat ion , 
the  WSA  with  quasi -crystal  line  approximation  iQCAl  yields  better  results. 

At  resonance  frequencies  we  note  that  P-YA  is  so  far  the  appropriate 
correlation  function  to  be  employed  [29-31], 


Formulation  of  the  Problem 

Consider  .V  identical,  finite  elastic  inclusions  that  are  randomly 

distributed  in  a  different  elastic  medium,  see  Fig.  1.  The  scatterers  are 

homogeneous  with  elastic  properties  given  by  Lame’s  constants  and  u,  and 

density  c ^ .  The  properties  of  the  outside  medium  (call  matrix)  are  given 

bv  X,  u  and  o.  In  Fig.  1,  0.  and  0.  refer  to  the  center  of  the  i-th  and  i-th 

i  j 

scatterers,  respectively  and  they  are  referred  to  the  origin  0  by  the 
spherical  polar  coordinates  (r^,  6^,  P  is  any  point  in  the  medium 

outside  the  scatterers  (the  matrix  medium) . 

A  time  harmonic  plane  wave  of  unit  amplitude  and  frequency  w  is 
incident  on  the  medium  such  that  the  direction  of  propagation  of  the  incident 
waves  is  along  the  c-axis,  which  may  be  written  in  terms  of  displacement 
field  vector  u^: 


-0 
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(r) 


i(k.=-  t)x 

s 

♦  e 


in 


where  and  k  are  the  compressions!  and  shear  wave  numbers  given  by 


k  =  m  /  c  ;  c  =  »  i  \  +  2  a  1  /  o 
P  P  ? 


k  =  uj/c  ;  c  - 
s  s  s 


and  t  is  the  time.  The  waves  incident  to  the  discrete  random  media  will 

->s  -*■ 

undergo  multiple  scattering.  Let  u^(r)  be  the  field  scattered  by  the  i-th 
scatterer.  The  incident  and  scattered  fields  satisfy  the  vector  Helmholtz 
equation.  The  problem  at  hand  reduces  to  computing  the  total  wave  field 
at  any  point  in  the  matrix  medium  and  hence  the  bulk  properties,  satisfying 
the  appropriate  boundary  condition  on  the  surface  of  the  scatterers  and 
radiation  conditions  at  infinity. 

The  total  field  at  any  point  in  the  matrix  medium  can  be  interpreted 
as  the  sum  of  the  incident  field  and  the  fields  scattered  by  all  the  scattere 
which  can  be  written  as 

u(r)  =  u  (r)  +  T  u*(pi)  ;  pi  =  r-ir  .  (4) 

i  =  1 


However,  the  field  that  excites  the  i-th  scatterer  is  the  incident  field  u 
plus  the  fields  scattered  from  all  other  scatterers  except  the  i-th.  The  ter: 
exciting  field  u  is  used  to  distinguish  between  the  field  actually  incident 
on  a  scatterer  and  the  external  incident  field  produced  by  a  source  at 
infinity.  Thus,  at  a  point  r  in  the  vicinity  of  the  i-th  scatterer,  write 
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where  'a'  is  a  typical  dimension  of  the  scatterer. 

The  exciting  and  scattered  fields  for  each  scatterer  can  be  expanded 
in  terms  of  rector  spherical  functions  with  respect  to  an  origin  at  the 
center  of  that  scatterer: 
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where  ^  _  ("  =  1.2,3)  are  the  vector  spherical  vector  basis  functions  [19]. 

Field  quantities  that  are  regular  at  the  origin  are  expanded  in  terms  of  the 

regular  (Re)  basis  set  (Re  )  obtained  by  replacing  the  Hankel  function 

of  the  first  kind,  h^,  in  the  above  equations  by  the  spherical  Bessel  functions 

j  of  the  first  kind.  In  Eq.  (7),  we  abbreviate  these  vector  basis  functions 

as  j  ,  =u  .  We  note  that  y,  is  for  the  longitudinal  part  while  i. 

’rino  tn  In  *  -ti 

^  id  u.  for  the  transverse  parts.  The  choice  of  the  basis  set  in  Eq.  1.7) 
on  r 

satisfies  the  radiation  condition  at  infinity  for  the  scattered  field, 
while  the  choice  in  Eq.  (6)  satisfies  the  regular  behavior  of  the  exciting 
field  in  the  region  a  <  jo. [  <  2a.  The  superscript  i  on  the  basis  functions 
refer  to  expansions  with  respect  to  CK  ,  and  and  B*  are  the  unknown 
exciting  and  scattered  field  coefficients.  We  also  expand  the  incident 
field  in  terms  of  vector  spherical  functions: 
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where  5  is  the  Kronecker  6.  For  the  sake  of  simplicity,  we  write  the 


incident  wave  field  in  terms  of  expansion  co-efficients  a  as  follows 

tn 


i  k  •  r . 

^  a  Re  +1  e  T  1 
tn  in 
tn 


(9) 


where  a  are  the  known  incident  field  coefficients, 
in 


The  unknown  coefficients  b1  can  be  related  to  B1  bv  means  of  anv 

in  in 


convenient  scattering  operator,  in  this  case  we  employ  the  T-matrix,  see 
Ref.  [52]. 


i1 

in 
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Substituting  Eqs.  (6),  (7)  and  (S)  in  (5),  we  obtain 
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Since  Che  field  quantities  are  expanded  with  respect  to  centers  of  each 
scatterer,  we  obtain  Eq.  (9)  with  basis  functions  with  respect  to  i-th 
and  j-th  centers.  In  order  to  express  them  with  respect  to  a  common 
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origin  (L  ,  we  employ  the  translation  and  addition  theorems  for  the  vector 
spherical  functions  [35]  which  may  be  written  in  a  compact  form  as  follows: 
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Employing  Eq.  (12)  in  (11)  and  using  the  orthogonality  of  the  vector 
spherical  basis  functions,  we  obtain  the  following  set  of  coupled  algebraic 
equations  for  the  exciting  field  coefficients  b* 
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x  'n' 
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(13) 


With  the  scattered  field  coefficients  expressed  in  terms  of  exciting 
field  coefficients  b^  and  the  T-matrix  as  given  by  (10),  Eq.  (13)  gives 
the  exciting  field  formulation  of  the  multiple  scattering.  If  we  multiply 
both  sides  of  Eq.  (15)  by  the  T-matrix,  then  we  obtain  the  scattered  field 
formulation  of  multiple  scattering  which  may  be  written  as 
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From  Eq.  (14),  it  can  be  seen  that  the  scattered  field  coefficients 
of  the  i-th  scatterer  explicitly  depend  on  the  position  and  orientation 
of  other  scatterers.  In  this  paper,  we  consider  a  random  distribution 
of  spherical  scatterers  and  the  case  when  N  -  °°  and  the  volume  occupied 
by  the  scatterers  V  -*■  »  such  that  N/V  =  n^  is  a  finite  number  density. 

For  such  distribution,  a  configurational  average  of  Eq .  (14)  can  be 
made  over  the  positions  of  all  scatterers  [28-52]  with  QCA  [24]  to  arrive 
at  an  equation  for  the  configurational  average  °f  the  scattered 

field  coefficients  with  one  scatterer  fixed: 
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where  p Cr^ | r ^ )  is  the  two  particle  joint  probability  density 
The  joint  probability  density  is  defined  as 


7  Sd^j-rj)  ;  |r.-r.  |  >  2a 


p(r . | r. ; 
K  J  i 


J  i  ~ 


r . -r.  >  2a 

J  i 


Equation  (16)  implies  that  the  particles  are  hard  (no- interpenet rat  ion )  and 
the  excluded  volume  is  a  sphere  of  radius  'a'  although  the  particles 
themselves  may  be  non  -  spherica 1 .  The  function  g(|r^-r^j)  is 
called  the  pair  correlation  function  and  depends  only  on  'r.-r. 1  due  to 
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translational  invariance  of  the  system  under  consideration.  The  pair 
correlation  function  for  an  ensemble  of  particles  depends  on  the  nature 
and  range  of  the  interparticle  forces.  The  average  of  several  measurements 
of  a  statistical  variable  that  characterises  an  ensemble  will  depend  on 
the  pair  correlation  function.  To  obtain  expressions  for  the  pair 
correlation  function,  one  needs  a  description  of  the  interparticle  forces. 
In  our  case  we  assume  that  the  scatterers  behave  like  effective  hard 
spheres  (where  the  radius  'a'  is  that  of  the  sphere  circumscribing  the 
scatterer) .  Percus  and  Yevick  [27]  have  obtained  an  approximate  integral 
equation  for  the  pair  correlation  function  of  a  classical  fluid  in 
equilibrium.  Wertheim  [34]  has  obtained  a  series  solution  of  the  integral 
equation  for  an  ensemble  of  hard  spheres.  The  statistics  of  the  fluid 
are  then  same  as  those  of  the  ensemble  of  discrete  hard  particles  that 
we  are  considering. 

Although  integral  expressions  for  the  correlation  functions  also 
result  in  a  heirarchy,  Percus  and  Yevick  have  truncated  the  heirarchy  by 
making  certain  approximations  that  result  in  a  self-consistent  relation 
between  the  pair  correlation  function  g(x)  and  the  direct  correlation 
function  C(x).  The  direct  correlation  function  may  be  interpreted  as 
the  correlation  function  resulting  from  an  'external  potential'  that 
produces  a  simultaneous  density  fluctuation  at  a  point  and  the  external 
potential  is  taken  to  be  the  potential  seen  by  a  particle  given  that 
there  is  a  particle  fixed  at  another  site.  Fisher  [35]  comments  that 
the  Percus-Yevick  approximation  is  a  strong  statement  of  the  extremely 
short  range  nature  of  the  direct  correlation  function.  The  integral 
equation  has  the  form 


...  -• 
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t(x)  =  1  *  n  j  x(x')dx'  -  n 

0  x<2a  0  x'<2a 


j  T (X ’ ) T (X -X ' ) dx ' 

< ; 

x-x ' | >2a 


where 


T(x)  =  g(x) 

x  > 

2a 

g(x)  =  0 

X  < 

2a 

T(X)  =  -C(X) 

X  < 

2a 

C(x)  =  0 

X  > 

2a 

(IS) 


Wertheim  [34]  has  solved  the  integral  equation  by  Laplace  transformation 
that  results  in  an  analytic  expression  for  C(x)  in  the  form 


C(x)  =  -(1-n)'4  [  ( l  +  2n) 2  -  6n(l*j  n)2x  *  n(H-2n)2  x-5/:]  ;  n  =  c/8 


where  'c'  is  the  effective  spherical  concentration  of  the  particles.  The 
Percus -Yevick  approximation  fails  as  the  concentration  approaches  the 
close  packing  factor  for  spheres  and  is  expected  to  be  good  for  c  <  0.3  or 
0.4. 

Equation  (19)  can  be  substituted  back  into  Eq .  (17)  to  yield  a 
series  solution  for  g(x)  in  the  form  [54] 


gU)  =  l  g  fx 
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where 


Sit)  = 


+  6n(l-n)t  +  18n“t-12n(l+2n) 


(22) 


and 


L(t)  =  1 2 n  [ ( l  +  n2) t  +  (l  +  2n)]  .  (23) 

Throop  and  Bearman  [36]  have  tabulated  g(x)  as  a  function  of  x  for  values 
of  n  =  c/8.  A  few  representative  plots  of  the  pair  correlation  function 
are  shown  in  Fig.  2. 

To  solve  the  integral  equations  given  by  (15),  we  consider  the 
inhomogeneous  medium  with  discrete  scatterers  as  a  homogeneous  continuum 
and  assume  that  the  average  coherent  wave  is  a  plane  wave  propagating 
with  an  effective  wave  number  K  in  the  same  direction  as  the  incident 
plane  wave,  h'e  can  thus  write 

i  K-r. 

<8*  >  =  X  e  1  (24) 

tn  in 


where  X  is  the  amplitude  of  the  coherent  wave, 
xn  r 

Substituting  Eq.  (24)  in  (15)  employing  the  joint  probability  function 
as  defined  before  and  the  divergence  theorem  to  convert  the  volume 
integral  in  (15)  to  surface  integrals  and  using  the  extinction  theorem 
which  cancels  the  incident  wave,  we  obtain  a  set  of  simultaneous  coupled 
homogeneous  equations  for  the  coefficients  X  given  by 


V 

L 

T"n’ 


n '  +n" 


i ' n '  q= j  n ' -n" J 


X  T 

x'n'  tn,T"n' 


CT'n' ,:"n" 


.  lL. 
(  k  ,  -  K ' 
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where  c  =  4m  a  n^/S  is  the  effective  spherical  concentration  of  the 
scatterers  per  unit  volume,  is  an  expression  containing  '.v'igner 
coefficients,  and 

I  (K, k  , c )  =  - ^ T  [2k  a  j  (2ka)  h'(2k  a) 

«  T  (kTa)2-(Ka)2  T  q  ^  T 

OD 

-2Ka  h  (2k  a)  j  (2ka)] -  24c  /  x2[g(x)-l]  h  (k_x)  j  (kx)  dx  (2b) 

4  T  4  x=l  4  4 

At  low  values  of  concentration  c,  g(x)  1,  see  Fig.  2,  and  hence  the 
integral  in  Eq.  (26)  is  negligible  which  results  in  a  system  of  uncorrelated 
hard  particle  statistics.  This  is  what  has  been  referred  to  as  the  'well 
stirred  approximation'  (1VSA)  and  yields  the  ’hole  correction  integral’  as 
outlined  by  Fikioris  and  Waterman  [37]  and  by  us  earlier.  If  g(x)  >  1, 
one  can  regard  the  Eq.  (26)  as  a  modified  ’hole  correction  integral'  which 
is  of  the  same  form  used  by  Twersky  [26], 

Equation  (25)  is  a  system  of  simultaneous  linear  homogeneous  equations 
for  the  unknown  amplitudes  X  .  For  nontrivial  solution,  we  require  that 
the  determinant  of  the  truncated  coefficient  matrix  vanishes,  which  yields 
an  equation  for  the  effective  wave  number  K  in  terms  of  k^  and  the  T-matrix 
of  the  scatterer.  This  is  the  dispersion  relation  for  the  scatterer 
filled  medium.  Equation  (25)  is  a  general  expression  valid  for  any 
arbitrary  shaped  scatterer,  since  the  T-matrix  is  the  only  factor  that 
contains  information  about  the  exact  shape  and  boundary  conditions  at  the 
scatterer.  Thus  the  formalism  presented  here  is  valid  for  all  the  three 
wave  fields.  The  effective  wave  number  K  obtained  in  the  analysis  is  a 


complex  quantity,  the  real  part  of  which  relates  to  the  phase  velocity, 
while  the  imaginary  part  relates  to  attenuation  of  coherent  waves  in  the 
medium . 


13 


Results  and  Conclusions 

In  the  Rayleigh  or  low  frequency  limit,  the  size  of  the  scatterers 
is  considered  to  be  small  when  compared  to  the  incident  wavelength.  It 
is  then  sufficient  to  take  only  the  lowest  order  coefficient  in  the 
expansion  of  the  fields.  In  this  limit,  the  elements  of  the  T-matrix  can 
be  obtained  in  closed  form  for  various  simple  shapes  (46).  It  can  be 
shown  that  at  low  frequencies,  only  *tq>  and  *T  ^  of  Eq.  (25) 
make  a  contribution.  After  some  manipulations  of  the  resulting 
3x3  determinant,  we  obtain  the  following  dispersion  relations  for 
elastic  spherical  inclusions  embedded  in  a  different  elastic  medium 
(matrix) : 


(U9c  El)(1  +  3c  Eq) 
1  -  15c  E,  £l  +  3c  Eq 


(22) 


k 

P 
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where 


.  3X+2u-(3\  +2u  ) 

p  =  —  - i - — 

0  3  4u+3X^+2p^ 


i  uOj-u)  24uJL(y1-u)  -  (,^  +  2^)  (191^  +  16*0 

24u  "(u^-'u I ~  (X~^7  (i9u1+'l6u) 


(29) 


_ 1 _ 

4y(u^-y)  +  3(A+2u)  (2'y^  +  3p) 

and  c  =  4^a~>n  /3  is  the  concentration  of  spheres,  and  K  and  are  the 
U  p  $ 

coherent  wave  numbers  for  longitudinal  and  shear  waves,  respectively,  in 
the  new  medium.  Similar  expressions  can  also  be  derived  for  spheroidal 
inclusions  using  the  T-matrix  obtained  in  Refs.  [32,38].  In  the  Rayleigh 
limit,  the  value  of  K  as  determined  by  the  above  dispersion  relations  is 
a  real  quantity  for  lossless  (elastic)  material  and  a  complex  quantity  for 
lossy  (viscoelastic)  material,  and  relates  to  phase  velocity  V  =  j/K. 

In  this  limit,  we  normally  study  the  dependence  of  phase  velocity  on 
concentration,  angle  of  incidence  and  aspect  ratio  o'7  *  scatterers.  The 
teneral  tendency  of  the  phase  velocity  is  to  in.  ost.  i  ghtly  (for 
inclusion)  and  decrease  slightly  (for  cracks  and  cavities)  as  concentration 
increases.  Thus,  the  phase  velocity  vs.  concentration  information  is  not 
very  useful  both  from  theoretical  3nd  experimental  point  of  view.  The 
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plots  of  absorption  and  coherent  attenuation  due  to  multiple  scattering 
vs.  frequency  for  various  concentrations  carry  more  information  which  may 
eventually  be  used  for  designing  absorbing  materials  [39], 


The  dispersion  relations  given  in  Eqs  .  (27)  and  (28)  may  also  be 
useful  in  obtaining  the  effective  shear  modulus  and  bulk  modulus  at  low 
frequencies.  Following  the  work  by  us  [20,22]  and  by  Bedeaux  and  Mazur  [2 
we  arrive  at  the  following  shear  and  bulk  moduli  (<y>  and  <B>)  of  an  elast 
material  containing  a  random  distribution  of  stress  free  bubbles  or 
cavities 


4y-3c  E,  (9\fl4u) 
4u+6c  E,  (5a+8u) 


(30) 


<B> 
B  ' 


3.\  +  2u  [  1  -6c  Eq] 
(3X+2u)  [ l+3c  Eq1 


(31) 


where  E^  and  E7  are  defined  in  Eq.  (29). 

To  study  the  response  at  resonant  and  higher  frequencies,  we  must 
consider  higher  powers  of  k^a,  and  this  implies  that  a  larger  number  of 
terms  (X  )  must  be  kept  in  the  expansion  of  the  average  field.  This  is 
best  done  numerically.  For  a  given  value  of  ka,  the  T-matrix  for  the 
scatterer  is  computed.  Next,  the  coefficient  matrix  M  corresponding 
to  X_n  (Eq.  (25)  is  formed.  The  complex  determinant  of  the 
coefficient  matrix  is  computed  using  standard  Gauss  elimination  techniques. 
For  a  given  k_a,  the  root  of  the  equation  det  M  =  0  is  searched  in  the 
complex  K  plane  (K^  +  ik,)  using  Muller's  method.  Good  initial  guesses 
were  provided  by  the  Rayleigh  limit  solutions  at  low  values  of  k^a  and 


16 


these  could  be  used  systematically  to  obtain  convergence  of  roots  at 
increasingly  higher  values  of  k^a.  The  real  part  K^  determines  the  phase 
velocity,  while  the  imaginary  part  K9  determines  the  coherent  wave 
attenuat ion . 


Here,  we  present  some  sample  numerical  calculation  of  spherical 
glass  inclusion  in  epoxy  matrix.  The  longitudinal  and  shear  wave 
velocities  of  the  glass  and  epoxy  matrix  are  taken  as  (c  )  =  5.28  mm/a  sec , 

c^  =  2.54  mm/ysec,  (cs^  =  3.24  mm/usec  and  c^  =  1.16  mm/usec,  respectively. 
IVe  consider  a  concentration  of  44.1}  to  reflect  a  high  concentration .  The 
coherent  wave  attenuation  vs.  frequency  (longitudinal  wave  number)  for 
this  configuration  is  shown  in  Fig.  5.  The  general  tendency  of  attenuation 
is  to  increase  at  lower  frequencies  and  shows  some  oscillation  as  shown. 
These  results  are  compared  with  some  experimental  observations  for  the 
same  composite  obtained  by  Kinra  (private  communication).  The  theoretical 
results  obtained  in  this  paper  compare  with  Kinra' s  experimental  results 
qualitatively  not  quantitatively.  The  reason  for  this  factor  difference 
must  be  explored  in  the  future.  The  oscillation  at  higher  frequencies, 
however,  indicate  that  the  scattering  is  mostly  in  the  forward  direction. 
Thus,  in  this  case  repeated  scattering  should  not  be  important,  since 
the  backscattered  wave  is  significantly  smaller  than  the  forward  scattered 
wave.  The  same  observation  may  be  noticed  even  for  electromagnetic  waves 
[28]  where  the  theoretical  results  obtained  by  our  theory  are  compared 
with  experimental  data.  (The  paper  [28]  is  enclosed  for  the  benefit  of 


the  reader. 


Since  the  phase  velocity  does  vary  very  slightly  as  a  function  of 
frequency,  the  bulk  properties  depend  totally  on  coherent  wave  at  ter.u.it :  on 


only.  Thus,  one  cun  compute  the  bulk  properties  which  cun  be  plotted 
in  the  complex  plane  (Cole-Cole  plot!  as  shown  in  our  paper  [50]  which 
is  also  enclosed. 

Recommendations  for  Future  ivork 

It  is  obvious  from  the  preceding  discussions  on  the  QCA  as  well  as 
the  numerical  results  that  the  two  major  improvements  required  are  tor 
the  QCA  as  well  as  the  pair  correlation  function,  so  that  good  results 
can  be  obtained  for  all  concentrations  even  at  long  and  intermediate 
values  of  the  wavelength.  In  a  review  article,  Lax  [40]  has  suggested 
that  in  the  quantum  mechanical  context,  the  QCA  could  be  improved  by 
using  modified  propagators  for  the  fields.  In  the  classical  context, 
this  implies  that  on  the  average,  single  particle  scattering  takes  place 
in  a  macroscopical ly  homogeneous  medium,  and,  in  this  respect,  this  idea 
is  the  same  as  the  coherent  potential  approximation  (CPA]  of  Solid  State 
Physics.  The  repeated  multiple  scattering  between  pairs  of  scatterers 
or  cluster  effects  can  be  improved  by  making  the  self  consistent 
approximation  (SCA)  in  addition  to  the  CPA. 

For  the  purpose  of  discussion  of  these  ideas  within  the  T-matrix 
formalism  given  earlier,  we  denote  by  u .  and  u^  the  tields  scattered  by 
and  exciting  the  j-th  scatterer,  respectively.  The  expansion  coefficien 
of  these  fields  are  denoted  by  B^  and  b'1 ,  respectively,  omitting  all 
subscripts . 

The  CPA  can  be  expressed  succinctly  as 

•  B'V  =  T(k'l  <b'V  ;5d 

l  J 


where  the  T- matrix  relating  the  exciting  and  scattered  field  coefficients 
is  evaluated  using  the  bulk  propagation  constant  K  for  the  embedding 
medium.  Thus  the  CPA  implies  that  the  field  scattered  by  a  single 
obstacle  in  the  presence  of  several  others  when  averaged  over  the  position 
of  all  seatterers  is  the  same  as  the  field  that  would  be  produced  by  a 
single  particle  embedded  in  a  macroscop ica 1 ly  homogeneous  medium  described 
by  the  propagation  constant  K.  The  incorporation  of  the  CPA  into  the 
previous  formalism  involves  changes  only  in  the  computations  and  a 
redefinition  of  the  T-matrix  in  Eq .  (10).  It  would  be  interesting  to  see 
the  change,  if  any,  in  the  numerical  computations  as  a  result  of  invoking 
the  CPA. 

The  idea  behind  the  'self  consistent  approximation'  (SCA)  is  somewhat 
more  suble.  From  the  discussion  in  the  section  on  the  QCA,  it  is  now  clear 
that  QCA -CPA  neglects  multiple  scattering  between  two  fixed  seatterers. 

The  SCA  as  defined  by  Schwart:  and  Ehrenreich  [41]  restores  this  by 
stating  that 


<BJ> . .  =  T(K)  <bj> . 
ij  J 


where  T(K)  is  the  T-matrix  of  scatterer  'j'  in  the  presence  of  scatterer 
'i'  in  the  effective  medium  with  propagation  constant  K.  Expressions  for 
T(K)  as  given  by  Varadan  and  Varadan  [42]  may  be  written  as 


R  t  r  L  .  /  2 )  T[l-3(-r.  .)Tj(r.  ^T] 


ri-H-r.  .  )TR(r.  .  )  ]  Rl-r.  ,/2) 

1  i  j  i  j  i  j 
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where  c  is  a  compact  notation  for  the  translation  matrices  B  and  C 
introduced  in  Eq.  [12).  The  R  matrix  is  simply  the  part  of  a  that  is 
regular  at  the  origin,  i.e.,  for  jr^.j  =  0.  All  matrices  in  Eq.  [341 
are  obtained  using  the  bulk  Propagation  Constant  for  the  host  medium. 

We  observe  that  TfK)  explicitly  depends  on  r^  the  distance  between 
' i '  and  'j1.  The  integration  procedure  will  no  longer  be  simple  as 
before  and  the  SCA  may  be  rather  difficult  to  enforce  in  computations, 
especially  if  more  realistic  models  are  chosen  for  the  pair  correlation 


tunct ion . 


Incc • porat ion  of  the  CPA  as  well  as  improved  models  of  the  pair 
correlation  function  into  our  computations  are  in  progress.  We  hope  that 
they  will  shed  some  light  on  the  sensitivity  of  multiple  scattering 
theories  to  approximations  like  QCA  and  SCA  as  a  function  of  frequency 
and  scatterer  concentration.  Needless  to  say  additional  experimental 
results  are  required  for  comparison  with  these  computations. 


I 
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Figure  1.  Random  distribution  of  inclusions  of  arbitrary  shape 
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Abstract 

Numerical  computations  of  the  effective  dielectric  constant  of  discrete  random  media 
are  presented  as  a  function  of  frequency.  Such  media  have  a  complex  dielectric 
constant  giving  rise  to  absorption  of  a  propagating  wave  both  due  to  geometric 
dispersion  or  multiple  scattering  as  well  as  absorption,  if  any,  due  to  the 
viscosity  of  the  particles  and  the  matrix  medium.  We  are  concerned  with  the 
absorption  due  to  multiple  scattering.  The  scattering  characteristics  of  the 
individual  particles  are  described  by  a  transition  or  T-matnx.  The  effects  of  two 
models  of  the  pair  correlation  function  which  arises  in  the  multiple  scattering 
analysis  are  considered,  he  conclude  that  the  well  stirred  approximation  iwSA)  is 
good  for  sparse  concentrations  and/or  high  frequencies  whereas  the  Percus-Yevick 
approximation  £?-YA)  is  preferred  for  higher  concentrations. 

Introduction 

The  study  of  the  frequency  dependence  of  the  effective  dielectric  constant  of 
statistically  inhomogeneous  media  is  important  for  practical  applications  such  as 
geopnystcal  exploration,  artificial  dielectrics  etc.  In  such  dielectrics  a 
propagating  electromagnetic  wave  undergoes  dispersion  and  absorption.  Some 
materials  are  naturally  absorptive  due  to  viscosity  whereas  inhomogeneous  media 
exhibit  absorption  due  to  geometric  dispersion  or  multiple  scattering. 

In  this  paper  the  effective,  complex  frequency  dependent  dielectric  constant  of  a 
discrete  random  medium  containing  a  distribution  of  ailgned  spheroidal  dielectric 
scatterers  in  free  space  is  calculated  for  different  concentrations  of  the  scatterers 
as  well  as  for  different  material  properties  of  the  scatterers.  we  use  a  multiple 
scattering  formalism  analogous  to  that  used  by  Twersky1  but  use  the  concept  of  a 
transition  matrix  or  T-matrix  to  characterice  the  scattering  from  a  single  obstacle. 
Ail  details  of  the  geometry  ana  material  properties  of  the  scatterer  are  contained 
in  the  T-matrix  leaving  the  general  formalism  indepenoent  of  the  type  of  scatterer. 
Spherical  statistics  are  used  even  though  the  scatterers  nav  be  non-spnencai .  Lax's 
quasi-crystalline  approximation  IQC.A)  is  used  to  truncate  the  beirarchv  of  (uations 
that  result  when  an  ensemble  average  is  performed  on  the  multiply  scattered  field. 


The  resulting  equations  tor  the  iverage  field  require  a  knowledge  or  the  pair 
correlation  function  of  the  dielectric  scatterers.  In  previous  work’’1'1'3.  we  assumed 
that  the  particles  did  not  penetrate  each  other  but  were  otherwise  uncorre latea. 
iillis’  has  called  this  the  well  stirred  approximation  , *SA).  However,  the  *SA  lead 
to  unpnysical  result,?  for  the  absorption  coefficient  of  the  average  medium  for 
scatterer  concentrations  c  >  0.123.  In  many  artificial  dielectrics,  the  scatterer 
concentration  is  often  greater  than  0.123.  In  this  Daper,  we  have  also  considered 
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the  Percus-Yevick  approximation  ;.P-YA)  to  the  pair  correlation  function.  Wertheim 

has  provided  a  semi-analytical  solution  of  the  resulting  integral  equation  for  a 
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system  at  hard  spheres.  Throop  and  Searman  have  provided  tabulated  values  or  the 
pair  correlation  function  for  different  values  of  the  concentration  as  a  function  of 
the  inter  particle  distance.  He  have  used  these  tabulated  values  in  the  numerical 
computations . 

Calculations  are  presented  for  a  system  of  polyethylene  spheres  and  spneroids  as  well 
as  ice  particles  for  0  <  c  <  0.26  for  several  values  of  the  non-dimensional 
wavenumber  ka  •  ^  ranging  from  0  to  3.0.  ('a'  is  a  characteristic  dimension  of  the 
scatterer).  Two  types  of  results  are  presented.  In  the  first  instance  the  validity 
of  the  WSA  and  P-YA  and  their  effect  on  the  absorption  coefficient  is  studied  as  a 
function  of  concentration  and  frequency.  Secondly,  the  complex  plane  locus  of  the 
effective  dielectric  constant  is  plotted  for  the  systems  considered.  For  artif-cial 
dielectrics  the  locus  deviates  dramatically  from  the  circular  arc  locus  commonlv 
noticed  for  ordinary  solids  and  liquids  that  exhibit  absorption  due  to  viscosity. 


*ave  propagation  in  a  discrete  random  medium 

Consider  N  identical  rotationally  symmetric  dielectric  scatterers  chat  are  aligned 
but  distributed  randomly  in  free  space  (see  Fig.  !:.  let  0  be  the  origin  of  a 
coordinate  system  located  outside  the  scatterers  whose  centers  are  denoted  '* 

0^,  0,,  0.  ...  0^.  Monochromatic  plane  electromagnetic  waves  of  frequency  j  propagate 
along  the  symmetry  axis  of  the  scatterers  which  is  taken  to  be  the  t-axis.  Since 
the  medium  is  isotropic  about  the  ;  -  axis  there  are  no  depoiantation  effects.  The 
time  dependence  of  the  incident  and  hence  the  fields  scattered  bv  the  maividua. 
scatterers  is  all  of  the  form  exp(-iut)  and  this  is  suppressed  in  the  equations  that 
follow. 

let  E0'.t)  be  the  electric  field  arising  from  the  incident  plane  wave  and  !*'r  '.he 
field  scattered  by  the  i-ch  scatterer.  The  total  field  at  a  point  r  outside  i..  the 
N  scatterers,  denoted  by  Eir)  is  given  by 
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I 
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The  field  incident  on  or  exciting  the  i-ch  scatterer  is  given  a 
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ana  the  vector  sortencal  functions  are  defined  as 
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In  Eqs.  (T)  and  (3)  ar.d  h^  are  the  spherical  Bessel  and  Hankel  functions  ana  the 
lre  *he  r,oraa'1-1:e^  spherical  harmonics  defined  with  real  angular  functions. 
To  make  the  notation  aore  compact  we  introduce  a  super  index  'n'  to  represent 
(timer  as  follows 


ou; 

Re 
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Re  "n 


We  observe  that  the  coefficients  of  exoansion  A*  and  F  associated  with  the 

•  n  n 

exciting  and  scattered  fields  depend  on  the  position  of  ail  N  scatterers.  Further, 
since  Eq.  ,;3)  is  satisfied,  we  can  relate  the  two  sets  of  coefficients  by  means  of 
the  T-aatnx  as  defined  by  Waterman1'0.  We  have 
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r*  a‘ 

nm  n 


The  T-matrix  depends  on  the  frequency  of  the  wave  exciting  the  scatterer  as  .ell  as 
its  geometry  and  material  properties. 

If  Eqs.  ‘,4) ,  (3)  and  v9)  are  substituted  m  Eq.  (t),  we  would  need  the  translation 
addition  theorems  for  the  vector  spherical  functions  m  order  that  we  aav  refer  all 
expansions  :n  Eq.  (2)  to  a  common  origin.  In  compact  form 
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wnere  r.  .  *  r.-r  is  the  vector  connecting  0.  to  0  and t  ,  is  the  translation  matrix 
L j  ii  i  i  nn ' 

for  the  vector  functions  and  R  ,  is  the  same  as  i  ,  with  the  spherical  Hankel 

nn '  nn '  r 

functions  in  j  reolaced  bv  spherical  Bessel  functions.  Detailed  expressions  for 
^  11 
the  matrices  are  given  by  Bostrom 

The  incident  electric  field  1°  can  be  expanded  with  respect  to  an  origin  at  0  as 


*o,-.  ik: 
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-mere  the  coefficients  are  known  ,see  for  example  Horse  ana  Fesncacn""' .  We 
oos erve  that  for  a  plane  wave  propagating  m  tne  indirection  tne  oniv  non-iero  values 
of  *  a-inc  are  and  4-»2,w  ie(l,*],  all  other  coefficients  being  :ero. 

'Jsmg  fas .  (4)  f  (S)  f  v9)-ill)  m  Eq.  '2',  using  the  orthogonality  of  tne  vector 
spherical  functions  we  obtain 
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Equation  -.12/  is  a  set  of  coupled  algebraic  equations  for  the  exciting  field 
coefficients  associated  with  each  scatterer.  If  the  number  of  scatterers  N  is 
finite  and  tne  position  of  ail  the  scatterers  is  known,  then  Ha.  .12,  can  be  solvec 
in  principle.  But  we  wish  to  consider  the  case  If-*,  V—*  sucn  that  N/Y  *  is  a 
finite  number  density.  Since  N  is  large,  we  are  oniy  interested  m  tne  conrigurationa 
average  of  Eq.  ;.12)  over  the  positions  of  11  particles 


The  conerent  field 

The  average  of  Eq.  '12'  over  the  position  of  ail  scatterers  the  average  exciting 
field)  is  the  sane  as  the  ensemble  average,  where  the  ensemble  is  composed  of 
different  possible  configurations  of  the  scatterers.  Equation  v 1 2 )  is  averaged  over 
the  position  of  all  particles  except  the  i-th.  3ut  the  right  hand  side  of  Eq.  12' 
explicitly  depends  on  the  position  of  the  ;-th  particle.  Hence  we  must  specify 
the  two  particle  joint  probability  density  ?■>.  rA  Further,  we  assume  that  all 
scatterers  are  identical,  so  that 


*e  note  that  the  average  exciting  field  with  one  scatterer  held  fixed  is  given  in  term 
or  the  average  with  two  scatterers  held  fixed,  leading  to  a  heirarcny  that  reauires 
knowledge  of  higher  order  probability  densities.  It  has  been  customary  to  truncate 
the  heirarchy  by  invoking  the  'quasi  crystalline  approximation'  i^QCA ;  first 
introduced  by  lax“.  According  to  this  approximation 


AA  ^  /v\\ 

V/  V/i 


specifically  the  QCA  neglects  multiple  scattering  between  pairs  of  scatterers. 
Improvements  to  the  QCA  have  been  suggested  by  Twersky*  and  in  previous  -.or*  :v  us'. 


The  joint  probability  density  is  defined  as 
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Euuacion  -  IS’i  implies  that  the  particles  ire  tara  no- interpenetration)  and  the 

excluded  voiume  is  a  spnere  or  radius  ’a*  although  the  particles  themselves  aav  -e 

non- spherical .  The  function  g(.  r.-t^  ;  is  caliea  the  oair  correlation  function  and 

depends  oniv  on  ri-ri  due  to  translational  invariance  of  the  system  inaer 
j  - 

consideration. 

we  assume  that  the  coherent  field  propagates  m  the  sane  direction  is  the  incident 
field  with  a  new,  effective  wavenumDer  <  that  is  Comdex  ana  freauencv  lepender.t . 
Hence 


where  A  is  the  amplitude  of  the  coherent  wave.  Thus  the  average  exciting  fie.d 
coefficient 


The  Kronecker  ielcas  in  £q.  Cl*'  indicate  that  only  tne  acimuthal  index  a*l 
contributes,  since  the  conerenc  wave  propagates  m  the  c-direction  and  those  in  the 
square  bracket  indicate  that  there  is  no  deco laritat ion. 

Equations  (14  -11')  are  substituted  in  Eq  (13'..  Since  the  T-matrix  of  i  rctationa 
symmetric  scatterer  is  block  diagonal  in  the  acimuthal  index  .see  Waterman ina 
the  coherent  field  propagates  in  the  :-direction,  the  sums  associatea  with  the 
atimuthal  indices  of  the  super  indices  n’  and  n"  m  Eq.  (13)  are  removed.  Further, 
as  shown  in  previous  work  bv  us3  as  well  as  Twersky',  the  extinction  theorem  can  oe 
used  to  cancel  the  incident  wave  term  in  Eo.  (13)  with  the  contribution  of  the 
integral  at  infinity.  Finally  Eq.  (13)  can  be  written  in  the  form 
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In  Eq.  (19)  c  *  ^5-  n  aJ  is  the  effective  spherical  concentration  of  the  particles 

[■?  . 

ll  12  lj"]  is  the  Wigner  J-;  symbol. 

BL|  »2  ®3j 

If  che  integral  in  Eq.  (19)  can  be  evaluated  for  suitable  models  of  the  pair 
correlation  function,  then  Eq.  (13)  is  a  set  of  coupled,  homogeneous,  algebraic 
equations  for  the  coherent  field  expansion  coefficients.  For  a  non-crivial  solution, 
the  determinant  of  the  coefficient  matrix  oust  vanish.  This  yields  the  required 
dispersion  equation  for  the  effective  or  average  medium.  In  general  che  system  of 
equations  can  be  solved  only  numerically  to  yield  the  effective  wave  number  11  as  a 
function  of  frequency  (k«Wc)  which  is  complex  (K  *  K^wiK^)  -  The  real  part  Kj 
yields  the  phase  velocity  in  the  medium  and  che  imaginary  part  IC,  leads  to  damping 
of  a  propagating  wave  due  to  geometric  dispersion  as  well  as  real  losses  if  any, 
associated  with  the  discrete  particles.  We  now  proceed  to  consider  the  evaluation 
of  the  integral  in  Eq.  (19). 


The  Percus-Ycvick  pair  correlation  function 

The  pair  correlation  function  for  an  ensemble  of  particles  depends  on  the  nature  and 
range  of  the  interparticle  forces.  The  average  of  several  measurements  of  a 
statistical  variable  that  characterizes  an  ensemble  will  depend  on  the  pair 
correlation  function.  As  we  have  seen,  the  coherent  or  average  electric  field  in  an 
ensemble  of  dielectric  scatterers  depends  on  the  pair  correlation  function  (Eqs.  .18)- 
(19)).  To  obtain  expressions  for  the  pair  correlation  function,  one  needs  a 
description  of  the  interparticle  forces.  In  our  case  we  assume  that  the  dielectric 
scatterers  behave  like  effective  hard  spheres  (where  the  radius  'a'  is  that  of  the 
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sphere  circumscribing  the  scatterer).  Percus  and  Yevick  have  obtained  an 
approximate  integral  equation  for  the  pair  correlation  function  of  a  classical  fluid 
in  equilibrium.  Wertheim8  has  obtained  a  series  solution  of  the  integral  equation 
for  an  ensemble  of  hard  spheres.  The  statistics  of  the  fluid  are  then  same  as  chose 
of  the  ensemble  of  discrete  hard  particles  that  we  are  considering. 

Although  Integral  expressions  for  the  correlation  functions  also  result  in  a 
heirarehy,  Percus  and  Yevick  have  truncated  the  heirarchy  by  making  certain 
approximations  that  result  in  a  self-consistent  relation  between  the  pair  correlation 
Amctlon  g(x)  and  the  direct  correlation  function  C(x).  The  direct  correlation 


function  may  be  interpreted  u  tfte  correlation  function  resulting  from  an  'external 
potential'  that  produces  a  simultaneous  density  fluctuation  at  a  point  and  the 
external  potential  is  taken  to  be  the  potential  seen  ov  a  particle  given  that  there 
is  a  particle  fixed  at  another  site.  Fisher12  comments  that  the  Percus-Yevick 
approximation  is  a  sfrong  statement  of  the  extremely  short  range  nature  of  the 
direct  correlation  function.  The  integral  equation  has  the  fora 
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where 


t (x)  ■  g(x)  ;  x  >  2a 
g(x)  «  0  ;  x  <  2a 

t(x)  ■  -C(x)  ;  x  <  2a 
C(x)  •  0  ;  x  >  2a 

a 

Nertheim  has  solved  the  integral  equation  by  Laplace  transformation  that  results  in 
an  analytic  expression  for  C(n)  in  the  form 


C(x)  •  -(1-n)*4  [ Cl-2n)2  -  6n(l^n)2x  *  n(l-2n)2  x3/2]  ;  n  •  c/S  (23) 


where  ’c'  is  ehe  effective  spherical  concentration  of  the  particles.  The  Percus- 
Yevick  approximation  fails  as  the  concentration  approaches  the  close  packing  factor 
for  spheres  and  is  expected  to  be  good  for  c  <  0.3  or  O.a. 

Equation  (23)  can  be  substituted  back  into  Eq.  (21)  to  yield  a  series  solution  for 
g(x)  in  the  form 
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Throop  and  Seaman  have  tabulated  g(x)  as  a  function  of  x  for  values  of  *  •  c/S.  A 
few  representative  plots  of  the  pair  correlation  function  are  shown  in  Fig.  2.  These 
tabulated  values  were  used  In  evaluating  the  integral  in  Eq.  (19). 


t 
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Figure  2.  The  Percus- Yevick  pair  Figure  3.  Coherent  attenuation  vs. 

correlation  function  g(x)  concentration  for  spherical 

ice  particles 


Coaparlson  of  WSA  and  P-YA 

The  homogeneous  system  of  algebraic  equation*  for  the  effective  exciting  field  were 
solved  numerically  for  two  different  node  Is  of  the  correlation  integral  I  appearing 
in  Eq.  (18).  In  eq.  (19)  if  the  second  tens  is  set  equal  to  zero,  we  just  have  a 
system  of  uncorrelated  hard  particles.  This  is  what  we  have  referred  to  as  the  well 
stirred  approximation  (WSA)6  earlier.  Computations  were  also  performed  by  numerically 
evaluating  the  integral  in  eq.  (19)  by  using  the  tabulated  values  of  the  Percus- 
Yevick  approximation  to  the  pair  correlation  functions  provided  by  T.iroop  and 
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Seeraan  . 

In  Fig.  3,  the  specific  damping  5^  ■  an  is  plotted  as  a  function  of  concentration 

for  a  random  distribution  of  numerical  ice  particles  (cr  •  3.168)  in  free  space  at 
ka  ■  O.SS.  The  WSA  agrees  with  the  P-YA  solution  only  up  to  concentrations 
C^O.OTS  and  then  there  is  a  marked  difference  and  the  WSA  fails  completely  at 
C  >  0.12S  leading  to  unphysical  resuics.  In  Fig.  a,  the  calculations  are  repeated 


Figure  4.  Coherent  attenuation  vs.  Figure  S.  Coherent  attenuation  vs. 

concentration  for  spherical  concentration  for 

ice  particles  polyethylene  spheres 


for  the  lue  systea  at  a  higher  value  of  ka  •  0.95.  Here  the  WSA  agrees  with  P-YA  up 
to  C  «  0.1  and  in  Fig.  S  siailar  calculations  were  perforaed  for  polyethylene  spheres 
(cr  »  2.26)  at  ka  »  4.62.  For  this  case  MSA  and  P-YA  results  agree  up  to  C  ■  3.15. 

Froa  these  results  it  would  appear  that  although  the  nSa  is  very  poor  at  higher 
scatterer  concentrations,  the  results  improve  dramatically  at  higher  values  of  ka, 
yielding  reasonably  good  results  for  higher  concentrations.  The  natural  explanation 
is  that  at  higher  values  of  ka,  aultiple  scattering  effects  between  pairs  of  particles 
becoae  smaller  and  thus  pair  correlation  effects  are  not  significant  and  the  OCA  also 
becomes  more  exact.  Sue  for  arbitrary  concentration  and  frequency  it  is  safer  to  use 
the  Percus-Yevick  approximation. 

The  effective  dielectric  constant 

Once  the  effective  complex  wavenumber  X  has  been  computed  by  solving  £q.  (IS) 
numerically,  we  can  proceed  further  and  evaluate  the  effective  dielectric  constant  of 
the  medium  which  is  also  complex  and  frequency  dependent.  In  the  usual  way,  the 
dielectric  constant  r*(w)  of  the  random  medium  is  defined  as 

n 

C*(u)  •  V  •€,(«)♦  iC,(«l) 

r  k*  1 

where  c,  and  t,  are  the  real  and  imaginary  parts  of  Che  dielectric  constant  and  the 
1  * 


«a*'.03 


Figure  6.  Complex  plane  locus  of  the 

effective  dielectric  constant 
for  a  system  of  polyethylene 
spheres 


Figure  7.  Complex  plane  locus  of  the 

effective  dielectric  constant 
for  a  system  of  spherical 
ice  particles 


subscript  on  c*  denotes  'relative  to  the  matrix  medium’.  The  real  part  c1  is  related 
to  the  refractive  index  and  phase  velocity  in  the  artifical  medium  and  the  imaginary 
part  accounts  for  the  dashing  in  the  medium.  In  real  materials,  the  damping  is 
intrinsic  to  the  system  and  is  due  to  macroscopic  viscosity  of  the  dielectric.  For 
the  artificial  or  effective  medium  under  consideration,  in  addition  to  natural 
losses  there  is  damping  due  to  geometric  dispersion  or  scattering. 
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Cole  and  Cole  have  given  a  convenient  representation  of  the  dispersion  and 
absorption  in  a  dielectric  by  means  of  an  Argand  diagram  or  a  plot  in  the  complex 
c-plane  of  Cj  versus  c^,  each  point  of  the  plot  being  characteristic  of  a  particular 
frequency.  For  many  types  of  loss  mechanisms,  the  locus  of  the  points  is  a  semi 
circle  with  its  center  on  the  real  axis  or  a  circular  arc.  In  Ref.  14,  the  complex 
dielectric  constant  of  several  liquids  and  solids  is  plotted  conforming  to  the 
circular  arc. 

In  the  present  case  the  complex  dielectric  constant  c(u)  corresponding  to  the 
effective  wavenumber  1C  of  the  effective  medium  is  studied  for  several  values  of  the 
frequency.  Overall  results  show  a  dramatic  deviation  from  the  circular  arc  locus. 
This  is  to  be  expected  since  the  medium  is  artificial. 


At  the  present  tine  there  are  no  experimental  results  available  to  verify  these 
calculations.  The  practical  applications  of  these  computations  are  many.  Such 
calculations  will  provide  reasonable  estimates  of  the  frequency  dependence 
dielectric  constant  as  a  function  of  particle  concentration,  site  and  shape  for 
inhomogeneous  media. 
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